NPL MCMCMH Software: Release 1 - Example random walk
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 random walk 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 random walk \n\n');
NPL MCMCMH Software: Release 1 National Physical Laboratory, UK April 2018 Example random walk
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)'; jumping = @(ad)jumprwg(ad, 2.4*L); % 2.4 is a scaling factor which makes the jumping distribution more % dispersed.
Specify quantiles
Q = [2.5; 50; 97.5];
Implement MCMC calculations
A0 = pars*ones(1, N) + L*randn(2, N); [S, aP, Rh, Ne, AA, IAA] = mcmcmh(target, jumping, M, N, M0, Q, A0);
Display percentage acceptance
fprintf('Mean percentage acceptance: %3g \n', mean(aP)); fprintf('\n')
Mean percentage acceptance: 22.7273
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.004218 Parameter 2: 1.004226
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: 10667 Parameter 2: 10649
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.179049 Parameter 2: 4.49284 Standard deviation Parameter 1: 0.377056 Parameter 2: 0.377069 2.5 percentile Parameter 1: -0.884899 Parameter 2: 3.72752 Median Parameter 1: -0.191559 Parameter 2: 4.50619 97.5 percentile Parameter 1: 0.585931 Parameter 2: 5.19878
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); aaa = (0:dela:aa(end)+5*dela)';
Generate figure
figure; bar(aa, nan) hold on if j == 1 plot(aaa,gampdf(aaa,m0a/2,2/m0a)) hold off xlabel('\alpha/1') ylabel('p(\alpha)/1') title('MCMC samples: \alpha') legend('Posterior','Prior') elseif j == 2 plot(aaa,gampdf(aaa,m0d/2,2/(m0d*s0d^2))) hold off xlabel('\delta/1') ylabel('p(\delta)/1') title('MCMC samples: \delta') legend('Posterior','Prior') 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