NPL MCM2MCMC Software: Release 1 - Gauge block example 3
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 3
The gauge block example involves the response model
where is the length of the gauge block at the reference temperature
,
is the temperature,
is the coefficient of thermal expansion and
.
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 3 \n\n');
NPL MCM2MCMC Software: Release 1 National Physical Laboratory, UK April 2018 Gauge block example 3
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 = 5500; N = 100; M0 = 500;
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; pS = 5; kappa = sqrt(2*pS+1)*sigma; Z1 = randn(2*pS,M,N); Z2 = randn(2*pS,M,N); C1 = sum(Z1.^2); C2 = sum(Z2.^2); B = C1./(C1+C2); Y = (z-kappa) + (2*kappa)*squeeze(B);
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.000064 Parameter 2: 1.000061 Parameter 3: 1.000002
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: 304595 Parameter 2: 310312 Parameter 3: 489785
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.361 102.764 Parameter 2: 20.0012 19.8664 Parameter 3: 0.100008 0.10002 Standard deviation MCM MCMC Parameter 1: 12.0941 12.2247 Parameter 2: 1.15485 1.1539 Parameter 3: 0.00577458 0.00577671 0 percentiles MCM MCMC Parameter 1: 77.7683 77.7683 Parameter 2: 18 18 Parameter 3: 0.09 0.09 2.5 percentiles MCM MCMC Parameter 1: 83.264 83.5463 Parameter 2: 18.101 18.0817 Parameter 3: 0.0905017 0.0905074 50 percentiles MCM MCMC Parameter 1: 99.9962 102.023 Parameter 2: 20.0011 19.7996 Parameter 3: 0.100012 0.100028 97.5 percentiles MCM MCMC Parameter 1: 124.313 124.844 Parameter 2: 21.9019 21.8805 Parameter 3: 0.109498 0.109503 100 percentiles MCM MCMC Parameter 1: 134.507 134.507 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