NPL MCM2MCMC Software: Release 1 - Exponential example 1
Contents
Introduction
This software is concerned with converting a sample from a Bayesian posterior distribution corresponding to a particular choice of prior distribution derived using the Monte Carlo method (MCM), to a Bayesian posterior corresponding to a preferred prior distribution. A Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm is used to undertake this conversion.
For a given measurement problem, the software performs the following tasks:
1. Generate an MCM sample from the distribution of the quantity of interest. This distribution can be thought of as a Bayesian posterior distribution corresponding to a particular choice of prior distribution.
2. Convert the MCM sample into an MCMC sample using a Metropolis-Hastings MCMC algorithm. The underlying or target distribution is the posterior distribution corresponding to a preferred prior distribution.
3. Analyse the convergence of the MCMC chains.
4. Extract and display summary information for the posterior distribution.
Exponential example 1
These calculations involve the response model
for the case
On this basis, is assigned the distribution
MCM/MCMC calculations are performed involving this response model assuming a uniform prior for . Details can be found in the user manual for the software.
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 MCM2MCMC Software: Release 1'); fprintf('\nNational Physical Laboratory, UK'); fprintf('\nApril 2018 \n'); fprintf('\nExponential example 1 \n\n');
NPL MCM2MCMC Software: Release 1 National Physical Laboratory, UK April 2018 Exponential example 1
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 MCM sample
Sample influence quantity
b = 2; sigmaB = 0.2; B = b + sigmaB*randn(M,N);
Sample response
z = 50*exp(-b); sigma = 0.2; Y = z + sigma*randn(M,N);
Create array to store MCM samples
A0 = zeros(M,N,2);
Generate MCM sample corresponding to a propagation of distributions
A0(:,:,1) = Y./exp(-B); A0(:,:,2) = B;
Convert MCM sample to MCMC sample
Record the derivative of response model function
D = exp(-B);
Implement MCMC calculations
[A,Rhat,Neff,IS,IA] = mcm2mcmc(A0,D,M0);
Analyse chain convergence
Calculate and display percentage acceptance
accept = round(100*sum(IA(:))/(M*N));
fprintf('Percentage acceptance: %3g \n\n',accept);
Percentage acceptance: 89
Display convergence indices
fprintf('Convergence indices \n') for i = 1:length(Rhat) fprintf('Parameter %1g: %12.6f \n',[i,Rhat(i)]) end fprintf('\n')
Convergence indices Parameter 1: 1.000260 Parameter 2: 1.000227
Display effective numbers of independent draws
fprintf('Effective number of independent draws \n') for i = 1:length(Neff) fprintf('Parameter %1g: %12g \n',[i,round(Neff(i))]) end fprintf('\n')
Effective number of independent draws Parameter 1: 65785 Parameter 2: 68769
Extract summary information for posterior distribution
Q = [0 2.5 50 97.5 100]'; [abar0,s0,aQ0,V0,AA0] = mcmcsum(A0,M0,Q); [abar,s,aQ,V,AA] = mcmcsum(A,M0,Q);
Display summary information for posterior distribution
fprintf('Summary information for posterior distribution \n') fprintf('Mean \n'); fprintf(' MCM MCMC \n'); for i = 1:length(abar0) fprintf('Parameter %1g: %12g % 12g \n',[i,abar0(i),abar(i)]) end fprintf('\n') fprintf('Standard deviation \n'); fprintf(' MCM MCMC \n'); for i = 1:length(abar0) fprintf('Parameter %1g: %12g %12g \n',[i,s0(i),s(i)]) end fprintf('\n') for j = 1:size(aQ0, 1) fprintf('%g percentiles \n',Q(j)); fprintf(' MCM MCMC \n'); for i = 1:length(abar0) fprintf('Parameter %1g: %12g %12g \n',[i,aQ0(j,i),aQ(j,i)]) end fprintf('\n') end
Summary information for posterior distribution Mean MCM MCMC Parameter 1: 51.044 53.1227 Parameter 2: 2.00046 2.04042 Standard deviation MCM MCMC Parameter 1: 10.4684 10.9087 Parameter 2: 0.200583 0.200711 0 percentiles MCM MCMC Parameter 1: 21.3298 22.863 Parameter 2: 1.15813 1.20847 2.5 percentiles MCM MCMC Parameter 1: 33.5937 35.0078 Parameter 2: 1.60838 1.64747 50 percentiles MCM MCMC Parameter 1: 49.9733 52.0134 Parameter 2: 2.00067 2.04016 97.5 percentiles MCM MCMC Parameter 1: 74.486 77.4853 Parameter 2: 2.39459 2.43461 100 percentiles MCM MCMC Parameter 1: 118.657 118.657 Parameter 2: 2.88971 2.88971
Generate histogram of samples
Assign number of bins
nB = 30;
for j = 1:2
Assign histogram for MCM samples
a0 = A0(:,:,j); a0 = a0(:); [na0,aa0] = hist(a0,nB); dela0 = aa0(2)-aa0(1); na0n = na0/(sum(na0)*dela0); yl1 = max(na0n);
Assign histogram for MCMC samples
[na,aa] = hist(AA(:,j),nB); dela = aa(2)-aa(1); nan = na/(sum(na)*dela); yl2 = max(nan); yl = max(yl1,yl2);
Generate figure
figure; bar(aa0,na0n,'w') hold on plot(aa,nan,'o') ylim([0,yl+0.2*yl]) if j == 1 xlabel('\alpha/1') ylabel('p(\alpha)/1') title('MCM and MCMC samples: \alpha') elseif j == 2 xlabel('\beta/1') ylabel('p(\beta)/1') title('MCM and MCMC samples: \beta') end legend('MCM','MCMC')


end
References
[1] A B Forbes, An MCMC algorithm based on GUM Supplement 1 for uncertainty evaluation, Measurement, 2012, DOI:10.1016/j.measurement.2012.01.018
[2] A B Forbes, A two-stage MCM/MCMC algorithm for uncertainty evaluation, in Advanced Mathematical and Computational Tools for Metrology IX, F Pavese et al., eds., World Scientific, Singpore, pp 159-170, 2012
[3] 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
MCM2MCMC Open-Source Software (Release 1.0): Software for measurement uncertainty evaluation
A B Forbes, K. Jagan and I M Smith, Data Science Group, National Physical Laboratory, UK
Copyright NPL Management Limited 2018