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

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

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

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');

В этом случае кажется, что объем выборки 1 000 более, чем достаточен дать хорошую точность для следующей средней оценки.

bHat = mean(trace)

bHat =

   -1.6931    4.2569

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

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

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