Создайте цепи MCMC для многомерного нормального распределения с помощью гамильтонового сэмплера Монте-Карло (HMC).

Определите количество параметров для выборки и их средства.

NumParams = 100;
means = randn(NumParams,1);
standevs = 0.1;

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

logpdf = @(theta)normalDistGrad(theta,means,standevs);

Выберите начальную точку дискретизатора. Создайте HMC-дискретизатор и настройте его параметры.

startpoint = randn(NumParams,1);
smp = hmcSampler(logpdf,startpoint);
smp = tuneSampler(smp);

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

NumChains  = 4;
chains     = cell(NumChains,1);
Burnin     = 500;
NumSamples = 2000;
for c = 1:NumChains
    if c == 1
        showOutput = 1;
    else
        showOutput = 0;
    end
    chains{c} = drawSamples(smp,'Burnin',Burnin,'NumSamples',NumSamples,...
        'Start',randn(size(startpoint)),'VerbosityLevel',showOutput,'NumPrint',500);
end

|==================================================================================|
|   ITER   |    LOG PDF    |  STEP SIZE  |  NUM STEPS  |  ACC RATIO  |  DIVERGENT  |
|==================================================================================|
|      500 |  8.450463e+01 |   4.776e-01 |           5 |   9.060e-01 |           0 |
|     1000 |  8.034444e+01 |   4.776e-01 |           9 |   8.810e-01 |           0 |
|     1500 |  9.156276e+01 |   4.776e-01 |           2 |   8.867e-01 |           0 |
|     2000 |  8.027782e+01 |   2.817e-02 |           6 |   8.890e-01 |           0 |
|     2500 |  9.892440e+01 |   4.648e-01 |           2 |   8.904e-01 |           0 |

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

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

for p = 1:3
    subplot(3,1,p);
    plot(chains{1}(:,p));
    ylabel(smp.VariableNames(p))
    axis tight
end
xlabel('Iteration')

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

function [lpdf,glpdf] = normalDistGrad(X,Mu,Sigma)
Z = (X - Mu)./Sigma;
lpdf = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
glpdf = -Z./Sigma;
end