diagnostics

Create MCMC chains using a Hamiltonian Monte Carlo (HMC) sampler and compute MCMC diagnostics.

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

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

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

startpoint = randn(3,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.

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

Compute MCMC diagnostics and display the results. Compare the true means in means with the column titled Mean in the MCMCdiagnostics table. The true posterior means are within a few Monte Carlo standard errors (MCSEs) of the estimated posterior means. The HMC sampler has accurately recovered the true means. Similarly, the estimated standard deviations in the column SD are very near the true standard deviations in 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

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