NPL MCMCMH Software: Release 1 - Example independence chain
Contents
Introduction
This software implements the Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm to generate samples from a user defined target distribution. The user can specify the jumping distribution used to draw proposal samples.
The software can draw multiple chains of samples simultaneously and can hence evaluate convergence statistics.
For a given measurement problem, the software performs the following tasks:
1. Specify the target and jumping distributions and their parameters.
2. Draw samples of desired size using the Metropolis-Hastings algorithm.
3. Analyse the convergence of the MCMC chains.
4. Extract and display summary information for the posterior distribution.
Example
Consider the model
,
It is assumed that is known.
and
are the quantities of interest and Gamma priors are assigned to them.
,
.
Calculations are performed for the case
A Gaussian independence chain jumping distribution is used to draw samples from the joint distribution of and
.
Samples of and
can be obtained by taking thier exponent.
Licence agreement
The software is provided with a software user licence agreement and the use of the software is subject to the terms laid out in that agreement. By running the software, the user accepts the terms of the agreement.
Main program
Identify the software to the user
fprintf('\nNPL MCMCMH Software: Release 1'); fprintf('\nNational Physical Laboratory, UK'); fprintf('\nApril 2018 \n'); fprintf('\nExample independence chain \n\n');
NPL MCMCMH Software: Release 1 National Physical Laboratory, UK April 2018 Example independence chain
Close all figures and clear all variables
close all clear format short
Set seed for random number generation
rng(1);
Assign chain length, number of chains and burn-in period
M = 1100; N = 100; M0 = 100;
Generate data
Length of data and prior information
m = 10; m0d = 10; s0d = 0.1; m0a = 5;
Generating and
delta = 150; alpha = 0.5;
Generating data with mean and standard deviation sy
sy = 1; y = alpha*delta + sy*randn(m, 1);
Specify Target distribution
target = @(ad)tar(ad, y, sy, m0d, s0d, m0a);
Evaluate the variance matrix of
from the Hessian Determining the MAP estimates of
and
a0 = [log(alpha); log(delta)]; map = @(ad)-target(ad); [pars, ~, ~, ~, ~, H] = fminunc(map, a0); V = inv(H);
Local minimum found. Optimization completed because the size of the gradient is less than the default value of the optimality tolerance. Computing finite-difference Hessian using objective function.
Specifying the jumping distribution
L = chol(V)'; jump = @()jumpicg(pars*ones(1, N), L);
Specify quantiles
Q = [2.5; 50; 97.5];
Implement MCMC calculations
A0 = jump(); [S, aP, Rh, Ne, AA, IAA] = mcmcmhic(target, jump, M, N, M0, Q, A0);
Display percentage acceptance
fprintf('Mean percentage acceptance: %3g \n', mean(aP)); fprintf('\n')
Mean percentage acceptance: 93.94
Display convergence indices
fprintf('Convergence indices \n') for i = 1:length(Rh) fprintf('Parameter %1g: %12.6f \n', [i, Rh(i)]) end fprintf('\n')
Convergence indices Parameter 1: 1.000139 Parameter 2: 1.000142
Display effective numbers of independent draws
fprintf('Effective number of independent draws \n') for i = 1:length(Ne) fprintf('Parameter %1g: %12g \n', [i, round(Ne(i))]) end fprintf('\n')
Effective number of independent draws Parameter 1: 78212 Parameter 2: 77947
Display summary information for posterior distribution
fprintf('Summary information for posterior distribution \n') fprintf('Mean \n'); for i = 1:2 fprintf('Parameter %1g: %12g \n', [i, S(1,i)]) end fprintf('\n') fprintf('Standard deviation \n'); for i = 1:2 fprintf('Parameter %1g: %12g \n', [i, S(2,i)]) end fprintf('\n') fprintf(' 2.5 percentile \n'); for i = 1:2 fprintf('Parameter %1g: %12g \n', [i, S(3,i)]) end fprintf('\n') fprintf(' Median \n'); for i = 1:2 fprintf('Parameter %1g: %12g \n', [i, S(4,i)]) end fprintf('\n') fprintf(' 97.5 percentile \n'); for i = 1:2 fprintf('Parameter %1g: %12g \n', [i, S(5,i)]) end fprintf('\n')
Summary information for posterior distribution Mean Parameter 1: -0.171375 Parameter 2: 4.48517 Standard deviation Parameter 1: 0.379967 Parameter 2: 0.37999 2.5 percentile Parameter 1: -0.889944 Parameter 2: 3.71759 Median Parameter 1: -0.180144 Parameter 2: 4.49413 97.5 percentile Parameter 1: 0.595797 Parameter 2: 5.20331
Generate histogram of samples
Assign number of bins
nB = 30;
for j = 1:2
Assign histogram for MCMC samples
Aa = exp(AA(:, :, j)); Aa = Aa(:); [na, aa] = hist(Aa, nB); dela = aa(2) - aa(1); nan = na/(sum(na)*dela);
Generate figure
figure; bar(aa, nan) if j == 1 xlabel('\alpha/1') ylabel('p(\alpha)/1') title('MCMC samples: \alpha') elseif j == 2 xlabel('\delta/1') ylabel('p(\delta)/1') title('MCMC samples: \delta') end


end
References
[1] Bayesian Data Analyis, A Gelman, J B Carlin, H S Stern and D B Rubin, Chapman & Hall/CRC, Boca Raton, 2004 (Section 11.6)
Acknowledgements
The work described here was supported by the National Measurement System Directorate of the UK Department for Department for Business, Energy and Industrial Strategy.
Program identifier
MCMCMH Open-Source Software (Release 1.0): Software for measurement uncertainty evaluation
A B Forbes and K. Jagan, Data Science Group, National Physical Laboratory, UK
Copyright NPL Management Limited 2018