Создайте цепи MCMC с помощью гамильтонова образца Monte Carlo (HMC) и вычислите диагностику MCMC.

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

means = [1;-2;2];
standevs = [1;2;0.5];
logpdf = @(theta)normalDistGrad(theta,means,standevs);

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

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

Извлеките образцы из задней плотности, используя несколько независимых цепей. Выберите различные, случайным образом распределенные начальные точки для каждой цепочки. Укажите количество выгорающих образцов, которые будут отброшены с начала цепочки Маркова, и количество образцов, которые будут сгенерированы после выгорания.

NumChains  = 4;
chains     = cell(NumChains,1);
Burnin     = 500;
NumSamples = 1000;
for c = 1:NumChains
    chains{c} = drawSamples(smp,'Burnin',Burnin,'NumSamples',NumSamples,...
        'Start',randn(size(startpoint)));
end

Вычислите диагностику MCMC и просмотрите результаты. Сравните истинные средства в means с заголовком столбца Mean в MCMCdiagnostics таблица. Истинные задние средства находятся в пределах нескольких стандартных ошибок Монте-Карло (MCSE) оцененных задних средств. Пробоотборник HMC точно восстановил истинные средства. Аналогичным образом, расчетные стандартные отклонения в столбце SD очень близко к истинным стандартным отклонениям в standev.

MCMCdiagnostics = diagnostics(smp,chains)

MCMCdiagnostics=3×8 table
     Name      Mean        MCSE        SD          Q5        Q95       ESS      RHat
    ______    _______    ________    _______    ________    ______    ______    ____

    {'x1'}     1.0038    0.016474    0.96164    -0.58601     2.563    3407.4     1  
    {'x2'}    -2.0435    0.034933      1.999     -5.3476    1.1851    3274.5     1  
    {'x3'}     1.9957    0.008209    0.49693      1.2036    2.8249    3664.5     1  

means

means = 3×1

     1
    -2
     2

standevs

standevs = 3×1

    1.0000
    2.0000
    0.5000

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