NPL MCM2MCMC Software: Release 1 - Exponential example 2
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 2
These calculations involve the response model
where and there is a prior constraint that
.
The calculations are performed for the case
MCM/MCMC calculations are performed involving this response model assuming a uniform prior for and
with
. 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 2 \n\n');
NPL MCM2MCMC Software: Release 1 National Physical Laboratory, UK April 2018 Exponential example 2
Close all figures and clear all variables
close all clear format short
Set seed for random number generation
rng(2);
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 + 0.2*randn(M,N);
Sample responses
z1 = 50*exp(-b); sigma = 0.2; Y1 = z1 + sigma*randn(M,N);
Sample responses
z2 = z1 + sigma; Y2 = z2 + sigma*randn(M,N);
Create array to store MCM samples
A0 = zeros(M,N,3);
Generate MCM sample corresponding to a propagation of distributions
A0(:,:,1) = Y1./exp(-B); A0(:,:,2) = (Y2-Y1)./exp(-B); A0(:,:,3) = B;
Convert MCM sample to MCMC sample
Record the absolute value of the determinant of the Jacobian matrix of the response model function
D = exp(-2*B);
Prior distribution imposing constraint
P00 = ones(M,N); L00 = A0(:,:,2) < 0; P00(L00) = 0;
Ensure that all chains start at a feasible point
L001 = P00(1,:) == 0;
if sum(L001) > 0
Find a chain with feasible starting point...
[maxP,imax ] = max(P00(1,:));
...and start all infeasible chains at the same point
A0(1,L001,1) = A0(1,imax,1); A0(1,L001,2) = A0(1,imax,2); A0(1,L001,3) = A0(1,imax,3); D(1,L001) = D(1,imax); P00(1,L001) = P00(1,imax);
end
Implement MCMC calculations
[A,Rhat,Neff,IS,IA] = mcm2mcmc(A0,D,M0,P00);
Analyse chain convergence
Calculate and display percentage acceptance
accept = round(100*sum(IA(:))/(M*N));
fprintf('Percentage acceptance: %3g \n\n',accept);
Percentage acceptance: 59
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.001341 Parameter 2: 1.001195 Parameter 3: 1.001214
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: 27225 Parameter 2: 29557 Parameter 3: 29230
Extract summary information for posterior distribution
Q = [0 2.5 50 97.5 100]'; [abar0,s0,aQ0,V0,AA0] = mcmcsum(A0,0,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: 50.9786 54.7496 Parameter 2: 1.50199 2.58642 Parameter 3: 1.99931 2.07941 Standard deviation MCM MCMC Parameter 1: 10.4339 11.1776 Parameter 2: 2.20778 1.82106 Parameter 3: 0.200458 0.200442 0 percentiles MCM MCMC Parameter 1: 20.7127 21.7048 Parameter 2: -11.5402 4.97924e-05 Parameter 3: 1.13089 1.17175 2.5 percentiles MCM MCMC Parameter 1: 33.5339 36.2271 Parameter 2: -2.70473 0.137033 Parameter 3: 1.60594 1.6891 50 percentiles MCM MCMC Parameter 1: 49.9025 53.6124 Parameter 2: 1.44076 2.26341 Parameter 3: 1.99902 2.07955 97.5 percentiles MCM MCMC Parameter 1: 74.4313 79.8809 Parameter 2: 6.0489 6.913 Parameter 3: 2.39334 2.47275 100 percentiles MCM MCMC Parameter 1: 126.24 126.24 Parameter 2: 13.9628 13.9628 Parameter 3: 2.91324 2.91324
Generate histogram of samples
Assign number of bins
nB = 30;
for j = 1:3
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(aa,nan,'w') hold on plot(aa0,na0n,'o') ylim([0,yl+0.2*yl]) if j == 1 xlabel('\alpha_1/1') ylabel('p(\alpha_1)/1') title('MCM and MCMC samples: \alpha_1') elseif j == 2 xlabel('\alpha_2/1') ylabel('p(\alpha_2)/1') title('MCM and MCMC samples: \alpha_2') elseif j == 3 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