Создайте цепочки 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.017202e+01 |   5.303e-01 |           4 |   8.820e-01 |           0 |
|     1000 |  7.756258e+01 |   5.303e-01 |           8 |   8.500e-01 |           0 |
|     1500 |  9.227175e+01 |   5.303e-01 |           2 |   8.547e-01 |           0 |
|     2000 |  8.073805e+01 |   3.128e-02 |           6 |   8.610e-01 |           0 |
|     2500 |  9.672952e+01 |   5.161e-01 |           2 |   8.600e-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