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

$$ \eta = \alpha[1 + \beta_2(\beta_1 - b_1) ], ~~~
   \alpha = \frac{\eta}{1 + \beta_2(\beta_1 - b_1)}, ~~~
   \zeta|\eta \sim {\textrm{Beta}}(\eta-\kappa, \eta+\kappa,
   p_{\rm S}, p_{\rm S}), $$

where $$ \alpha $$ is the length of the gauge block at the reference temperature $$ b_1 $$, $$ \beta_1 $$ is the temperature, $$ \beta_2 $$ is the coefficient of thermal expansion and $$ \kappa = \sigma (2p_{\rm S}+1)^{1/2} $$.

The calculations are performed for the case

$$ b_1 = 20, ~~~ \sigma = 2, ~~~
   \zeta = 100, ~~~ p_{\rm S} = 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 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 $$ \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 Beta}(\eta-\kappa, \eta+\kappa, p_{\rm S},    p_{\rm S}) $$

  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