Линейная регрессионая модель

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

Моделируйте среднее Гауссово распределение как функцию предикторов и параметров модели как

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

,

.

Чтобы использовать 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 ®, который возвращает логарифм продукта априорного и вероятностного, и градиент этого логарифма. The logPosterior функция определяется в конце этого примера. Затем вызовите функцию с аргументами, чтобы задать logpdf входной параметр в hmcSampler функция.

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

Создайте HMC Sampler

Определите начальную точку, с которой будет начата выборка, и затем вызовите 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 (maximum-a-posteriori) апостериорной плотности. Вы можете начать выборку с любой точки, но часто более эффективно оценить точку 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, StepSize, и 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

Визуализация выборок

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

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

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

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

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

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

The 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