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.
Choose a starting point. Create the HMC sampler and tune its parameters.
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.
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=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
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.