NPL MCM2MCMC Software: Release 1 - Gauge block 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.
Gauge block example 2
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 2 \n\n');
NPL MCM2MCMC Software: Release 1 National Physical Laboratory, UK April 2018 Gauge block example 2
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; nu = 5; Z = randn(nu,M,N); C = sum(Z.^2); T = randn(M,N).*sqrt(nu*ones(M,N)./squeeze(C)); Y = z + sigma*squeeze(T);
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.000065 Parameter 2: 1.000075 Parameter 3: 1.000000
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: 302475 Parameter 2: 286029 Parameter 3: 500000
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.358 102.769 Parameter 2: 20.0012 19.8655 Parameter 3: 0.100008 0.100018 Standard deviation MCM MCMC Parameter 1: 12.2116 12.3392 Parameter 2: 1.15485 1.15363 Parameter 3: 0.00577458 0.00577382 0 percentiles MCM MCMC Parameter 1: 40.6442 40.6442 Parameter 2: 18 18 Parameter 3: 0.09 0.09 2.5 percentiles MCM MCMC Parameter 1: 82.9839 83.2996 Parameter 2: 18.101 18.0819 Parameter 3: 0.0905017 0.0905063 50 percentiles MCM MCMC Parameter 1: 99.9946 102.026 Parameter 2: 20.0011 19.7983 Parameter 3: 0.100012 0.100034 97.5 percentiles MCM MCMC Parameter 1: 124.594 125.113 Parameter 2: 21.9019 21.8801 Parameter 3: 0.109498 0.109499 100 percentiles MCM MCMC Parameter 1: 169.899 169.899 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