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