NPL MCM2MCMC Software: Release 1 - Exponential 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.

Exponential example 1

These calculations involve the response model

$$ \eta = \alpha {\rm e}^{-\beta}, ~~~ \zeta|\eta \sim {\rm N}(\eta,\sigma^2), ~~~ \beta \sim {\rm N}(b,\sigma_{\rm B}^2), $$

for the case

$$ b = 2, ~~~ \sigma = 0.2, ~~~ \sigma_{\rm B} = 0.2, ~~~ \zeta = 50 {\rm e}^{-b}. $$

On this basis, $$ \eta|\zeta $$ is assigned the distribution $$ \eta|\zeta \sim {\rm N}(\zeta,\sigma^2). $$

MCM/MCMC calculations are performed involving this response model 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('\nExponential example 1 \n\n');

NPL MCM2MCMC Software: Release 1
National Physical Laboratory, UK
April 2018 

Exponential example 1

Close all figures and clear all variables

  close all
  clear
  format short

Set seed for random number generation

  rng(1);

Assign chain length, number of chains and burn-in period

  M  = 1100;
  N  = 100;
  M0 = 100;

Generate MCM sample

Sample influence quantity $$ \beta \sim {\rm N}(b, \sigma_{\rm B}^2) $$

  b = 2;
  sigmaB = 0.2;
  B = b + sigmaB*randn(M,N);

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

  z = 50*exp(-b);
  sigma = 0.2;
  Y = z + sigma*randn(M,N);

Create array to store MCM samples

  A0 = zeros(M,N,2);

Generate MCM sample corresponding to a propagation of distributions

  A0(:,:,1) = Y./exp(-B);
  A0(:,:,2) = B;

Convert MCM sample to MCMC sample

Record the derivative of response model function

  D = exp(-B);

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:  89

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.000260 
Parameter 2:     1.000227

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:        65785 
Parameter 2:        68769

Extract summary information for posterior distribution

  Q = [0 2.5 50 97.5 100]';
  [abar0,s0,aQ0,V0,AA0] = mcmcsum(A0,M0,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:       51.044      53.1227 
Parameter 2:      2.00046      2.04042 

Standard deviation 
                  MCM          MCMC 
Parameter 1:      10.4684      10.9087 
Parameter 2:     0.200583     0.200711 

0 percentiles 
                  MCM          MCMC 
Parameter 1:      21.3298       22.863 
Parameter 2:      1.15813      1.20847 

2.5 percentiles 
                  MCM          MCMC 
Parameter 1:      33.5937      35.0078 
Parameter 2:      1.60838      1.64747 

50 percentiles 
                  MCM          MCMC 
Parameter 1:      49.9733      52.0134 
Parameter 2:      2.00067      2.04016 

97.5 percentiles 
                  MCM          MCMC 
Parameter 1:       74.486      77.4853 
Parameter 2:      2.39459      2.43461 

100 percentiles 
                  MCM          MCMC 
Parameter 1:      118.657      118.657 
Parameter 2:      2.88971      2.88971

Generate histogram of samples

Assign number of bins

  nB = 30;

  for j = 1:2

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/1')
      ylabel('p(\alpha)/1')
      title('MCM and MCMC samples: \alpha')
    elseif j == 2
      xlabel('\beta/1')
      ylabel('p(\beta)/1')
      title('MCM and MCMC samples: \beta')
    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


Published with MATLAB® R2017b