Байесовский анализ для модели логистической регрессии

В этом примере показано, как сделать байесовские выводы для логистической регрессионой модели с помощью slicesample.

Статистические выводы обычно основаны на максимальной оценке правдоподобия (MLE). MLE выбирает параметры, которые максимизируют вероятность данных, и интуитивно более привлекательны. В MLE параметры считаются неизвестными, но фиксированными и оцениваются с некоторым доверием. В байесовской статистике неопределенность относительно неизвестных параметров количественно определяется с помощью вероятности, так что неизвестные параметры рассматриваются как случайные переменные.

Байесовский вывод

Байесовский вывод - это процесс анализа статистических моделей с включением предшествующих знаний о параметрах модели или модели. Корнем такого вывода является теорема Байеса:

$P(\mathrm{parameters|data}) = \frac
 {P(\mathrm{data|parameters}) \times P(\mathrm{parameters})}
 {P(\mathrm{data})}
 \propto \mathrm {likelihood} \times \mathrm{prior}$

Например, предположим, что у нас есть нормальные наблюдения

$X|\theta \sim N(\theta, \sigma^2)$

где sigma известна, а предшествующее распределение для theta определяется следующим образом:

$\theta \sim N(\mu, \tau^2)$

В этой формуле также известны mu и tau, иногда известные как гиперпараметры. Если мы наблюдаем n выборки X, мы можем получить апостериорное распределение для theta как

$\theta|X \sim N\left(\frac{\tau^2}{\sigma^2/n + \tau^2} \bar X +
 \frac{\sigma^2/n}{{\sigma^2/n + \tau^2}} \mu,
 \frac{(\sigma^2/n)\times \tau^2}{\sigma^2/n +
 \tau^2}\right)$

На следующем графике показаны априорные, вероятностные и апостериорные значения theta.

rng(0,'twister');

n = 20;
sigma = 50;
x = normrnd(10,sigma,n,1);
mu = 30;
tau = 20;
theta = linspace(-40, 100, 500);
y1 = normpdf(mean(x),theta,sigma/sqrt(n));
y2 = normpdf(theta,mu,tau);
postMean = tau^2*mean(x)/(tau^2+sigma^2/n) + sigma^2*mu/n/(tau^2+sigma^2/n);
postSD = sqrt(tau^2*sigma^2/n/(tau^2+sigma^2/n));
y3 = normpdf(theta, postMean,postSD);
plot(theta,y1,'-', theta,y2,'--', theta,y3,'-.')
legend('Likelihood','Prior','Posterior')
xlabel('\theta')

Данные автомобильного эксперимента

В некоторых простых задачах, таких как предыдущий пример вывода нормального среднего, легко вычислить апостериорное распределение в закрытой форме. Но в общих задачах, связанных с несопряженными априорными распределениями, апостериорные распределения трудно или невозможно вычислить аналитически. В качестве примера рассмотрим логистическую регрессию. Этот пример включает эксперимент, чтобы помочь смоделировать долю автомобилей с различными весами, которые не проходят тест пробега. Данные включают наблюдения веса, количества протестированных автомобилей и количества неудачных. Мы будем работать с преобразованной версией весов, чтобы уменьшить корреляцию в наших оценках параметров регрессии.

% A set of car weights
weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
weight = (weight-2800)/1000;     % recenter and rescale
% The number of cars tested at each weight
total = [48 42 31 34 31 21 23 23 21 16 17 21]';
% The number of cars that have poor mpg performances at each weight
poor = [1 2 0 3 8 8 14 17 19 15 17 21]';

Модель логистической регрессии

Логистическая регрессия, частный случай обобщенной линейной модели, подходит для этих данных, поскольку переменная отклика является биномиальной. Модель логистической регрессии может быть записана как:

$P(\mathrm{failure}) = \frac{e^{Xb}}{1+e^{Xb}}$

где X - матрица проекта, а b - вектор, содержащий параметры модели. В MATLAB ® мы можем записать это уравнение как:

logitp = @(b,x) exp(b(1)+b(2).*x)./(1+exp(b(1)+b(2).*x));

Если у вас есть некоторые предшествующие знания или некоторые неинформативные априорные значения доступны, можно задать априорные распределения вероятностей для параметров модели. Например, в этом примере мы используем нормальные априорные значения для точки пересечения b1 и наклонные b2, т.е.

prior1 = @(b1) normpdf(b1,0,20);    % prior for intercept
prior2 = @(b2) normpdf(b2,0,20);    % prior for slope

По теореме Байеса совместное апостериорное распределение параметров модели пропорционально продукту вероятностей и априоров.

post = @(b) prod(binopdf(poor,total,logitp(b,weight))) ...  % likelihood
            * prior1(b(1)) * prior2(b(2));                  % priors

Обратите внимание, что нормализующая константа для апостериорной функции в этой модели аналитически неразрешима. Однако, даже не зная нормирующей константы, вы можете визуализировать апостериорное распределение, если знаете приблизительную область значений параметров модели.

b1 = linspace(-2.5, -1, 50);
b2 = linspace(3, 5.5, 50);
simpost = zeros(50,50);
for i = 1:length(b1)
    for j = 1:length(b2)
        simpost(i,j) = post([b1(i), b2(j)]);
    end;
end;
mesh(b2,b1,simpost)
xlabel('Slope')
ylabel('Intercept')
zlabel('Posterior density')
view(-110,30)

Этот апостериор вытянут вдоль диагонали в пространстве параметров, что указывает на то, что после того, как мы посмотрим на данные, мы считаем, что параметры коррелируют. Это интересно, так как до того, как мы собрали любые данные, мы предполагали, что они независимы. Корреляция происходит от объединения нашего предыдущего распределения с функцией правдоподобия.

Отбор проб среза

Методы Монте-Карло часто используются в байесовском анализе данных, чтобы суммировать апостериорное распределение. Идея в том, что, даже если вы не можете вычислить апостериорное распределение аналитически, можно сгенерировать случайную выборку из распределения и использовать эти случайные значения для оценки апостериорного распределения или производной статистики, такой как апостериорное среднее, медиана, стандартное отклонение и т. Д. Выборка среза является алгоритмом, разработанным для выборки из распределения с произвольной функцией плотности, известной только до константы пропорциональности - именно то, что необходимо для выборки из сложного апостериорного распределения, константа нормализации которого неизвестна. Алгоритм не генерирует независимые выборки, а скорее марковскую последовательность, стационарное распределение которой является целевым распределением. Таким образом, срез sampler является алгоритмом Маркова Цепи Монте-Карло (MCMC). Однако он отличается от других хорошо известных алгоритмов MCMC, потому что необходимо указать только масштабированный апостериор - не нужно никаких предложений или маргинальных распределений.

В этом примере показов, как использовать выборку среза в рамках байесовского анализа модели логистической регрессии пробега, включая генерацию случайной выборки из апостериорного распределения для параметров модели, анализ выхода выборки и выводы о параметрах модели. Первый шаг - сгенерировать случайную выборку.

initial = [1 1];
nsamples = 1000;
trace = slicesample(initial,nsamples,'pdf',post,'width',[20 2]);

Анализ выходных данных сэмплера

После получения случайной выборки из селектора среза важно исследовать такие проблемы, как сходимость и смешивание, чтобы определить, может ли выборка разумно рассматриваться как набор случайных реализаций из целевого апостериорного распределения. Просмотр маргинальных графиков является самым простым способом исследования выхода.

subplot(2,1,1)
plot(trace(:,1))
ylabel('Intercept');
subplot(2,1,2)
plot(trace(:,2))
ylabel('Slope');
xlabel('Sample Number');

Из этих графиков очевидно, что эффекты начальных значений параметра занимают некоторое время, чтобы исчезнуть (возможно, пятьдесят или около того выборок), прежде чем процесс начнет выглядеть стационарным.

Также полезно при проверке на сходимость использовать скользящее окно для вычисления статистики, такой как среднее значение выборки, медиана или стандартное отклонение для выборки. Это создает более плавный график, чем следы необработанных выборок, и может облегчить идентификацию и понимание любой нестационарности.

movavg = filter( (1/50)*ones(50,1), 1, trace);
subplot(2,1,1)
plot(movavg(:,1))
xlabel('Number of samples')
ylabel('Means of the intercept');
subplot(2,1,2)
plot(movavg(:,2))
xlabel('Number of samples')
ylabel('Means of the slope');

Поскольку это скользящие средние значения в окне 50 итераций, первые 50 значений не сопоставимы с остальной частью графика. Однако остальная часть каждого графика, по-видимому, подтверждает, что апостериорное средство параметра сходилось к стационарности после 100 или около того итераций. Также очевидно, что эти два параметра коррелируют друг с другом в соответствии с более ранним графиком апостериорной плотности.

Поскольку период осаждения представляет выборки, которые не могут быть разумно обработаны как случайные реализации из целевого распределения, вероятно, целесообразно не использовать первые 50 или более значений в начале выходных данных селектора среза. Можно просто удалить эти строки выхода, однако также можно задать период «burn-in». Это удобно, когда подходящая длина ожога уже известна, возможно, по предыдущим запускам.

trace = slicesample(initial,nsamples,'pdf',post, ...
                                     'width',[20 2],'burnin',50);
subplot(2,1,1)
plot(trace(:,1))
ylabel('Intercept');
subplot(2,1,2)
plot(trace(:,2))
ylabel('Slope');

Эти графики трассировки, по-видимому, не показывают никакой нестационарности, что указывает на то, что период горения сделал свою работу.

Однако существует второй аспект трассировочных графиков, которые также должны быть исследованы. Хотя трассировка для точки пересечения выглядит как высокая частота шум, трассировка для наклона, по-видимому, имеет более низкую частотную составляющую, что указывает на автокорреляцию между значениями при смежных итерациях. Мы все еще могли вычислить среднее значение от этого автокоррелированной выборки, но часто удобно уменьшить требования к хранению, удаляя избыточность в выборке. Если бы это исключило автокорреляцию, это также позволило бы нам рассматривать это как выборку независимых значений. Для примера можно утончить образец, сохранив только каждое 10-е значение.

trace = slicesample(initial,nsamples,'pdf',post,'width',[20 2], ...
                                                'burnin',50,'thin',10);

Чтобы проверить эффект этого утончения, полезно оценить функции автокорреляции выборки из трасс и использовать их, чтобы проверить, быстро ли смешиваются выборки.

F    =  fft(detrend(trace,'constant'));
F    =  F .* conj(F);
ACF  =  ifft(F);
ACF  =  ACF(1:21,:);                          % Retain lags up to 20.
ACF  =  real([ACF(1:21,1) ./ ACF(1,1) ...
             ACF(1:21,2) ./ ACF(1,2)]);       % Normalize.
bounds  =  sqrt(1/nsamples) * [2 ; -2];       % 95% CI for iid normal

labs = {'Sample ACF for intercept','Sample ACF for slope' };
for i = 1:2
    subplot(2,1,i)
    lineHandles  =  stem(0:20, ACF(:,i) , 'filled' , 'r-o');
    lineHandles.MarkerSize = 4;
    grid('on')
    xlabel('Lag')
    ylabel(labs{i})
    hold on
    plot([0.5 0.5 ; 20 20] , [bounds([1 1]) bounds([2 2])] , '-b');
    plot([0 20] , [0 0] , '-k');
    hold off
    a  =  axis;
    axis([a(1:3) 1]);
end

Значения автокорреляции при первой задержке значительны для параметра точки пересечения и тем более для параметра наклона. Мы могли бы повторить выборку, используя больший параметр утончения в порядок, чтобы уменьшить корреляцию дальше. Однако в целях этого примера мы продолжим использовать текущую выборку.

Вывод для параметров модели

Как ожидалось, гистограмма выборки имитирует график апостериорной плотности.

subplot(1,1,1)
hist3(trace,[25,25]);
xlabel('Intercept')
ylabel('Slope')
zlabel('Posterior density')
view(-110,30)

Можно использовать гистограмму или оценку плотности сглаживания ядра, чтобы суммировать предельные свойства распределения апостериорных выборок.

subplot(2,1,1)
hist(trace(:,1))
xlabel('Intercept');
subplot(2,1,2)
ksdensity(trace(:,2))
xlabel('Slope');

Можно также вычислить описательную статистику, такую как апостериорное среднее или процентилы из случайных выборок. Чтобы определить, является ли размер выборки достаточно большим для достижения желаемой точности, полезно контролировать желаемую статистику трасс как функцию от количества выборок.

csum = cumsum(trace);
subplot(2,1,1)
plot(csum(:,1)'./(1:nsamples))
xlabel('Number of samples')
ylabel('Means of the intercept');
subplot(2,1,2)
plot(csum(:,2)'./(1:nsamples))
xlabel('Number of samples')
ylabel('Means of the slope');

В этом случае представляется, что размер выборки 1000 более чем достаточен, чтобы дать хорошую точность для апостериорной средней оценки.

bHat = mean(trace)

bHat =

   -1.6931    4.2569

Сводные данные

Statistics and Machine Learning Toolbox™ предлагает множество функций, которые позволяют легко задавать вероятности и априорные значения. Они могут быть объединены, чтобы вывести апостериорное распределение. The slicesample функция позволяет вам провести байесовский анализ в MATLAB с помощью симуляции Марковской цепи Монте-Карло. Он может использоваться даже в задачах с апостериорными распределениями, которые трудна выборка от использования стандартных генераторов случайных чисел.

Документация