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

В этом примере показано, как сделать байесовские выводы для модели логистической регрессии с помощью 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)$

где сигма известен и предшествующее распределение для тета

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

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

$\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)$

Следующий график показывает предшествующий, правдоподобный и задний для тета.

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)

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

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

Методы Монте-Карло часто используются в байесовском анализе данных для суммирования заднего распределения. Идея заключается в том, что, даже если вы не можете вычислить апостериорное распределение аналитически, вы можете сгенерировать случайную выборку из распределения и использовать эти случайные значения для оценки апостериорного распределения или производной статистики, такой как апостериорное среднее, медиана, стандартное отклонение и т.д. Выборка срезов - это алгоритм, предназначенный для выборки из распределения с произвольной функцией плотности, известной только до константы пропорциональности - именно то, что необходимо для выборки из сложного заднего распределения, константа нормализации которого неизвестна. Алгоритм генерирует не независимые выборки, а марковскую последовательность, стационарное распределение которой является целевым распределением. Таким образом, пробоотборник среза является алгоритмом Маркова цепочки Монте-Карло (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 или около того значений в начале выходных данных дискретизатора среза. Эти строки вывода можно просто удалить, однако можно также указать период «записи». Это удобно, когда подходящая длина горения уже известна, возможно, из предыдущих прогонов.

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

Резюме

Toolbox™ статистики и машинного обучения предлагает множество функций, которые позволяют легко определять возможности и приоры. Они могут быть объединены для получения заднего распределения. slicesample функция позволяет выполнять байесовский анализ в MATLAB с помощью моделирования Monte Carlo цепи Маркова. Его можно использовать даже при проблемах с задними распределениями, которые трудно отсчитать от использования стандартных генераторов случайных чисел.

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