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

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

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

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

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

,

.

Чтобы использовать 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, 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.6200

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

Чертите выборки от следующей плотности, с помощью нескольких независимых цепочек. Выберите различные начальные точки для цепочек, случайным образом распределенных вокруг предполагаемой точки 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.483179e+02 |   2.955e-01 |          10 |   9.200e-01 |           0 |
|      600 | -1.490021e+02 |   2.484e-02 |           3 |   9.117e-01 |           0 |
|      900 | -1.509802e+02 |   2.534e-01 |           4 |   9.144e-01 |           0 |
|     1200 | -1.496713e+02 |   1.317e-01 |          14 |   9.058e-01 |           0 |
|     1500 | -1.486536e+02 |   2.955e-01 |          10 |   9.113e-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.3855     0.005599    0.27813     1.9187     2.8472    2467.6    1.0006
    'x2'       2.5494    0.0063242    0.33267     1.9995      3.097    2767.1    1.0002
    'x3'     -0.44917    0.0066201    0.34165    -1.0259    0.10024    2663.4    1.0001
    'x4'    -0.060541    0.0028155    0.14199    -0.2831    0.18424    2543.4    1.0002


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