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

$$ \eta = \alpha[1 + \beta_2(\beta_1 - b_1) ], ~~~
   \alpha = \frac{\eta}{1 + \beta_2(\beta_1 - b_1)}, ~~~
   \zeta|\eta \sim {\rm N}(\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. $$

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 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 $$ \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 N}(\zeta, \sigma^2) $$

  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