Create MCMC chains for a multivariate normal distribution using a Hamiltonian Monte Carlo (HMC) sampler.

Define the number of parameters to sample and their means.

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

First, save a function normalDistGrad on the MATLAB® path that returns the multivariate normal log probability density and its gradient (normalDistGrad is defined at the end of this example). Then, call the function with arguments to define the logpdf input argument to the hmcSampler function.

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

Choose a starting point of the sampler. Create the HMC sampler and tune its parameters.

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

Draw samples from the posterior density, using a few independent chains. Choose different, randomly distributed starting points for each chain. Specify the number of burn-in samples to discard from the beginning of the Markov chain and the number of samples to generate after the burn-in. Set the 'VerbosityLevel' to print details during sampling for the first chain.

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 |

After obtaining a random sample, investigate issues such as convergence and mixing to determine whether the samples represent a reasonable set of random realizations from the target distribution. To examine the output, plot the trace plots of the samples for the first few variables using the first chain.

A number of burn-in samples have been removed to reduce the effect of the sampling starting point. Furthermore, the trace plots look like high-frequency noise, without any visible long-range correlation between the samples. This indicates that the chain is well mixed.

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

The normalDistGrad function returns the logarithm of the multivariate normal probability density with means in Mu and standard deviations in Sigma, specified as scalars or columns vectors the same length as the startpoint. The second output argument is the corresponding gradient.

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