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

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

В этой формуле также известны mu и tau, иногда известном как гиперпараметры. Если мы наблюдаем 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');

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

bHat = mean(trace)

bHat =

   -1.6931    4.2569

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

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

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