Создайте цепи 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')

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