Модель линейной регрессии

Анализируйте модель линейной регрессии с прерыванием, линейные коэффициенты (вектор-столбец), и шумовое отклонение распределения данных как свободные параметры. Примите, что каждая точка данных имеет независимое Распределение Гаусса:

Смоделируйте среднее значение Распределения Гаусса как функция предикторов и параметров модели как

В Байесовом анализе также необходимо присвоить предшествующие распределения всем свободным параметрам. Присвойте независимое Гауссово уголовное прошлое на прерывании и линейных коэффициентах:

,

.

Чтобы использовать HMC, все переменные выборки должны быть неограничены, означая, что следующая плотность и ее градиент должны быть четко определены для всех действительных значений параметров. Если у вас есть параметр, который ограничивается к интервалу, то необходимо преобразовать этот параметр в неограниченный. Чтобы сохранить вероятность, необходимо умножить предшествующее распределение на соответствующий якобиевский фактор. Кроме того, примите этот фактор во внимание при вычислении градиента следующего.

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

.

Запишите логарифм следующей плотности свободных параметров как

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

Создайте набор данных

Задайте истинные значения параметров для прерывания, линейные коэффициенты Beta, и шумовое стандартное отклонение. Знание истинных значений параметров позволяет соответствовать выходу сэмплера HMC. Только первый предиктор влияет на ответ.

NumPredictors = 2;

trueIntercept = 2;
trueBeta = [3;0];
trueNoiseSigma = 1;

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

NumData = 100;
rng('default') %For reproducibility
X = rand(NumData,NumPredictors);
mu = X*trueBeta + trueIntercept;
y = normrnd(mu,trueNoiseSigma);

Задайте плотность апостериорной вероятности

Выберите средние значения и стандартные отклонения Гауссова уголовного прошлого.

InterceptPriorMean = 0;
InterceptPriorSigma = 10;
BetaPriorMean = 0;
BetaPriorSigma = 10;
LogNoiseVarianceMean = 0;
LogNoiseVarianceSigma = 2;

Сохраните функциональный logPosterior на пути MATLAB®, который возвращает логарифм продукта предшествующего и вероятности и градиента этого логарифма. logPosterior функция задана в конце этого примера. Затем вызовите функцию с аргументами, чтобы задать logpdf входной параметр к hmcSampler функция.

logpdf = @(Parameters)logPosterior(Parameters,X,y, ...
    InterceptPriorMean,InterceptPriorSigma, ...
    BetaPriorMean,BetaPriorSigma, ...
    LogNoiseVarianceMean,LogNoiseVarianceSigma);

Создайте сэмплер HMC

Задайте начальную точку, чтобы начать произвести от, и затем вызывать hmcSampler функция, чтобы создать гамильтонов сэмплер как HamiltonianSampler объект. Отобразите свойства сэмплера.

Intercept = randn;
Beta = randn(NumPredictors,1);
LogNoiseVariance = randn;
startpoint = [Intercept;Beta;LogNoiseVariance];
smp = hmcSampler(logpdf,startpoint,'NumSteps',50);

smp

smp = 

  HamiltonianSampler with properties:

                  StepSize: 0.1000
                  NumSteps: 50
                MassVector: [4x1 double]
              JitterMethod: 'jitter-both'
      StepSizeTuningMethod: 'dual-averaging'
    MassVectorTuningMethod: 'iterative-sampling'
                    LogPDF: [function_handle]
             VariableNames: {4x1 cell}
                StartPoint: [4x1 double]

Оцените точку MAP

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

[MAPpars,fitInfo] = estimateMAP(smp,'VerbosityLevel',0);
MAPIntercept = MAPpars(1)
MAPBeta = MAPpars(2:end-1)
MAPLogNoiseVariance = MAPpars(end)

MAPIntercept =

    2.3857


MAPBeta =

    2.5495
   -0.4508


MAPLogNoiseVariance =

   -0.1007

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

plot(fitInfo.Iteration,fitInfo.Objective,'ro-');
xlabel('Iteration');
ylabel('Negative log density');

Настройте сэмплер

Важно выбрать хорошие значения для параметров сэмплера, чтобы получить эффективную выборку. Лучший способ найти хорошие значения состоит в том, чтобы автоматически настроить MassVectorНеродной размер, и NumSteps параметры с помощью tuneSampler метод. Используйте метод для:

1. Настройте MassVector из сэмплера.

2. Настройте StepSize и NumSteps в течение фиксированной продолжительности симуляции, чтобы достигнуть определенного приемного отношения. Целевое приемное отношение по умолчанию 0,65 хорошо в большинстве случаев.

Начните настраиваться в предполагаемой точке MAP для более эффективной настройки.

[smp,tuneinfo] = tuneSampler(smp,'Start',MAPpars);

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

figure;
plot(tuneinfo.StepSizeTuningInfo.StepSizeProfile);
xlabel('Iteration');
ylabel('Step size');

accratio = tuneinfo.StepSizeTuningInfo.AcceptanceRatio

accratio =

    0.6400

Чертите выборки

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

Установите 'VerbosityLevel' значение, чтобы распечатать детали во время выборки для первой цепи.

NumChains = 4;
chains = cell(NumChains,1);
Burnin = 500;
NumSamples = 1000;
for c = 1:NumChains
    if (c == 1)
        level = 1;
    else
        level = 0;
    end
    chains{c} = drawSamples(smp,'Start',MAPpars + randn(size(MAPpars)), ...
        'Burnin',Burnin,'NumSamples',NumSamples, ...
        'VerbosityLevel',level,'NumPrint',300);
end

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      300 | -1.484164e+02 |   2.532e-01 |          12 |   9.500e-01 |           0 |
|      600 | -1.492436e+02 |   2.128e-02 |           4 |   9.450e-01 |           0 |
|      900 | -1.509753e+02 |   2.171e-01 |           5 |   9.444e-01 |           0 |
|     1200 | -1.493455e+02 |   1.128e-01 |          16 |   9.358e-01 |           0 |
|     1500 | -1.489602e+02 |   2.532e-01 |          12 |   9.373e-01 |           0 |

Исследуйте диагностику сходимости

Используйте diagnostics метод, чтобы вычислить стандартную диагностику MCMC. Для каждого параметра выборки метод использует все цепи, чтобы вычислить эти статистические данные:

Отобразите таблицу диагностики и истинные значения параметров выборки, заданных в начале примера. Поскольку предшествующее распределение неинформативно для этого набора данных, истинные значения между (или рядом) 5-е и 95-е квантили.

diags = diagnostics(smp,chains)
truePars = [trueIntercept;trueBeta;log(trueNoiseSigma^2)]

diags =

  4x8 table

     Name       Mean         MCSE         SD          Q5         Q95       ESS       RHat 
    ______    _________    _________    _______    ________    _______    ______    ______

    {'x1'}       2.3805    0.0053354    0.28028      1.8966     2.8335    2759.7    1.0005
    {'x2'}       2.5544    0.0062264    0.33478      1.9989     3.1026    2890.9    1.0004
    {'x3'}     -0.44516    0.0064468     0.3415     -1.0247    0.10198      2806         1
    {'x4'}    -0.062826    0.0028183    0.14112    -0.28345    0.17997    2507.5    1.0004


truePars =

     2
     3
     0
     0

Визуализируйте выборки

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

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

figure;
subplot(2,2,1);
plot(chains{1}(:,1));
title(sprintf('Intercept, Chain 1'));
for p = 2:1+NumPredictors
    subplot(2,2,p);
    plot(chains{1}(:,p));
    title(sprintf('Beta(%d), Chain 1',p-1));
end
subplot(2,2,4);
plot(chains{1}(:,end));
title(sprintf('LogNoiseVariance, Chain 1'));

Объедините цепи в одну матрицу и создайте графики рассеивания и гистограммы, чтобы визуализировать 1D и 2D крайние апостериорные распределения.

concatenatedSamples = vertcat(chains{:});
figure;
plotmatrix(concatenatedSamples);
title('All Chains Combined');

Функции для вычисления апостериорного распределения

logPosterior функция возвращает логарифм продукта нормальной вероятности и нормального предшествующего для линейной модели. Входной параметр Parameter имеет формат [Intercept;Beta;LogNoiseVariance]X и Y содержите значения предикторов и ответа, соответственно.

normalPrior функция возвращает логарифм многомерной нормальной плотности вероятности со средними значениями Mu и стандартные отклонения SigmaВ виде скаляров или столбцов векторизовал ту же длину как P. Вторым выходным аргументом является соответствующий градиент.

function [logpdf, gradlogpdf] = logPosterior(Parameters,X,Y, ...
    InterceptPriorMean,InterceptPriorSigma, ...
    BetaPriorMean,BetaPriorSigma, ...
    LogNoiseVarianceMean,LogNoiseVarianceSigma)


% Unpack the parameter vector
Intercept        = Parameters(1);
Beta             = Parameters(2:end-1);
LogNoiseVariance = Parameters(end);
% Compute the log likelihood and its gradient
Sigma                   = sqrt(exp(LogNoiseVariance));
Mu                      = X*Beta + Intercept;
Z                       = (Y - Mu)/Sigma;
loglik                  = sum(-log(Sigma) - .5*log(2*pi) - .5*Z.^2);
gradIntercept1          = sum(Z/Sigma);
gradBeta1               = X'*Z/Sigma;
gradLogNoiseVariance1	= sum(-.5 + .5*(Z.^2));
% Compute log priors and gradients
[LPIntercept, gradIntercept2]           = normalPrior(Intercept,InterceptPriorMean,InterceptPriorSigma);
[LPBeta, gradBeta2]                     = normalPrior(Beta,BetaPriorMean,BetaPriorSigma);
[LPLogNoiseVar, gradLogNoiseVariance2]  = normalPrior(LogNoiseVariance,LogNoiseVarianceMean,LogNoiseVarianceSigma);
logprior                                = LPIntercept + LPBeta + LPLogNoiseVar;
% Return the log posterior and its gradient
logpdf               = loglik + logprior;
gradIntercept        = gradIntercept1 + gradIntercept2;
gradBeta             = gradBeta1 + gradBeta2;
gradLogNoiseVariance = gradLogNoiseVariance1 + gradLogNoiseVariance2;
gradlogpdf           = [gradIntercept;gradBeta;gradLogNoiseVariance];
end

function [logpdf,gradlogpdf] = normalPrior(P,Mu,Sigma)
Z          = (P - Mu)./Sigma;
logpdf     = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
gradlogpdf = -Z./Sigma;
end