NPL MCM2MCMC Software: Release 1 - Gauge block 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.
Gauge block example 1
The gauge block example involves the response model
where is the length of the gauge block at the reference temperature
,
is the temperature, and
is the coefficient of thermal expansion.
The calculations are performed for the case
MCM/MCMC calculations are performed involving this response model for simulating measurements of a gauge block 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('\nGauge block example 1 \n\n');
NPL MCM2MCMC Software: Release 1 National Physical Laboratory, UK April 2018 Gauge block example 1
Close all figures and clear all variables
close all clear format short
Set seed for random number generation
rng(3);
Assign chain length, number of chains and burn-in period
M = 1100; N = 100; M0 = 100;
Generate MCM sample
Sample temperature
b1 = 20; B1 = b1 + 2*(2*rand(M,N)-1);
Sample coefficient of thermal expansion
b2 = 0.1; B2 = b2 + 0.01*(2*rand(M,N)-1);
Sample response
z = 100; sigma = 2; Y = z + 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) = Y./(1 + B2.*(B1-b1)); A0(:,:,2) = B1; A0(:,:,3) = B2;
Convert MCM sample to MCMC sample
Record the derivative of response model function
D = abs(1 + B2.*(B1-b1));
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: 93
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.000139 Parameter 2: 1.000116 Parameter 3: 1.000065
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: 78276 Parameter 2: 81243 Parameter 3: 88533
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: 101.333 102.742 Parameter 2: 20.0046 19.8691 Parameter 3: 0.0999784 0.0999807 Standard deviation MCM MCMC Parameter 1: 12.0897 12.2173 Parameter 2: 1.15601 1.15481 Parameter 3: 0.00577961 0.00578233 0 percentiles MCM MCMC Parameter 1: 76.4781 77.8602 Parameter 2: 18.0001 18.0001 Parameter 3: 0.0900003 0.0900003 2.5 percentiles MCM MCMC Parameter 1: 83.2705 83.5541 Parameter 2: 18.102 18.084 Parameter 3: 0.0905127 0.0905162 50 percentiles MCM MCMC Parameter 1: 99.9574 102.055 Parameter 2: 20.0042 19.7962 Parameter 3: 0.0999586 0.0999513 97.5 percentiles MCM MCMC Parameter 1: 124.3 124.815 Parameter 2: 21.9014 21.8821 Parameter 3: 0.109519 0.109517 100 percentiles MCM MCMC Parameter 1: 134.506 134.506 Parameter 2: 22 22 Parameter 3: 0.11 0.11
Generate histograms 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(aa0,na0n,'w') hold on plot(aa,nan,'o') ylim([0,yl+0.2*yl]) if j == 1 xlabel({'\alpha/mm'}) ylabel({'p(\alpha)/mm^{-1}'}) title('MCM and MCMC samples: gauge block length \alpha') elseif j == 2 xlabel({'\beta_1/\circC'}) ylabel({'p(\beta_1)/\circC^{-1}'}) title('MCM and MCMC samples: temperature \beta_1') elseif j == 3 xlabel({'\beta_2/\circC^{-1}'}) ylabel({'p(\beta_2)/\circC'}) title('MCM and MCMC samples: coefficient of thermal expansion \beta_2') 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