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

$$ \eta = \alpha[1 + \beta_2(\beta_1 - b_1) ], ~~~
   \alpha = \frac{\eta}{1 + \beta_2(\beta_1 - b_1)}, ~~~
   \zeta|\eta \sim {\rm t}_{\nu}(\eta,\sigma^2), $$

where $$ \alpha $$ is the length of the gauge block at the reference temperature $$ b_1 $$, $$ \beta_1 $$ is the temperature, and $$ \beta_2 $$ is the coefficient of thermal expansion.

The calculations are performed for the case

$$ b_1 = 20, ~~~ \sigma = 2, ~~~
   \zeta = 100, ~~~ \nu = 5. $$

MCM/MCMC calculations are performed involving this response model for simulating measurements of a gauge block assuming a uniform prior for $$ \alpha $$. 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 $$ \beta_1 \sim {\rm R}(18,22) $$

  b1 = 20;
  B1 = b1 + 2*(2*rand(M,N)-1);

Sample coefficient of thermal expansion $$ \beta_2 \sim {\rm R}(0.09,0.11) $$

  b2 = 0.1;
  B2 = b2 + 0.01*(2*rand(M,N)-1);

Sample response $$ \eta|\zeta \sim {\rm t}_{\nu}(\zeta, \sigma^2) $$

  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