NPL MCMCMH Software: Release 1 - Example random walk

Contents

Introduction

This software implements the Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm to generate samples from a user defined target distribution. The user can specify the jumping distribution used to draw proposal samples.

The software can draw multiple chains of samples simultaneously and can hence evaluate convergence statistics.

For a given measurement problem, the software performs the following tasks:

1. Specify the target and jumping distributions and their parameters.

2. Draw samples of desired size using the Metropolis-Hastings algorithm.

3. Analyse the convergence of the MCMC chains.

4. Extract and display summary information for the posterior distribution.

Example

Consider the model

$$ y_i = \alpha \delta + \epsilon_i$,

$$\epsilon_i \sim N(0,\sigma_y^2), i = 1, ..., m.$$

It is assumed that $\sigma_y$ is known. $\alpha$ and $\delta$ are the quantities of interest and Gamma priors are assigned to them.

$$ \alpha \sim G(m_{0,a}/2, m_{0,a}/2)$,

$$ \delta \sim G(m_{0,d}/2, m_{0,d}\sigma_{0,d}^2/2)$.

Calculations are performed for the case

$$m_{0,a} = 5,~~~ m_{0,d} = 10,~~~, \sigma_{0,d} = 0.1,~~~ \sigma_y =
1.$$

A Gaussian random walk jumping distribution is used to draw samples from the joint distribution of $\log(\alpha)$ and $\log(\delta)$.

Samples of $\alpha$ and $\delta$ can be obtained by taking thier exponent.

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 MCMCMH Software: Release 1');
  fprintf('\nNational Physical Laboratory, UK');
  fprintf('\nApril 2018 \n');
  fprintf('\nExample random walk \n\n');
NPL MCMCMH Software: Release 1
National Physical Laboratory, UK
April 2018 

Example random walk 

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 data

Length of data and prior information

  m = 10;
  m0d = 10;
  s0d = 0.1;
  m0a = 5;

Generating $\delta$ and $\alpha$

  delta = 150;
  alpha = 0.5;

Generating data with mean $\delta \alpha$ and standard deviation sy

  sy = 1;
  y = alpha*delta + sy*randn(m, 1);

Specify Target distribution

  target = @(ad)tar(ad, y, sy, m0d, s0d, m0a);

Evaluate the variance matrix of $\log(\alpha)$ $\log(\delta)$ from the Hessian Determining the MAP estimates of $\log(\alpha)$ and $\log(\delta)$

  a0 = [log(alpha); log(delta)];
  map = @(ad)-target(ad);
  [pars, ~, ~, ~, ~, H] = fminunc(map, a0);
  V = inv(H);
Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.




Computing finite-difference Hessian using objective function.

Specifying the jumping distribution

  L = chol(V)';
  jumping = @(ad)jumprwg(ad, 2.4*L);
% 2.4 is a scaling factor which makes the jumping distribution more
% dispersed.

Specify quantiles

  Q = [2.5; 50; 97.5];

Implement MCMC calculations

  A0 = pars*ones(1, N) + L*randn(2, N);
  [S, aP, Rh, Ne, AA, IAA] = mcmcmh(target, jumping, M, N, M0, Q, A0);

Display percentage acceptance

  fprintf('Mean percentage acceptance: %3g \n', mean(aP));
  fprintf('\n')
Mean percentage acceptance: 22.7273 

Display convergence indices

  fprintf('Convergence indices \n')
  for i = 1:length(Rh)
      fprintf('Parameter %1g: %12.6f \n', [i, Rh(i)])
  end
  fprintf('\n')
Convergence indices 
Parameter 1:     1.004218 
Parameter 2:     1.004226 

Display effective numbers of independent draws

  fprintf('Effective number of independent draws \n')
  for i = 1:length(Ne)
      fprintf('Parameter %1g: %12g \n', [i, round(Ne(i))])
  end
  fprintf('\n')
Effective number of independent draws 
Parameter 1:        10667 
Parameter 2:        10649 

Display summary information for posterior distribution

  fprintf('Summary information for posterior distribution \n')
  fprintf('Mean \n');
  for i = 1:2
      fprintf('Parameter %1g: %12g \n', [i, S(1,i)])
  end
  fprintf('\n')

  fprintf('Standard deviation \n');
  for i = 1:2
     fprintf('Parameter %1g: %12g \n', [i, S(2,i)])
  end
  fprintf('\n')


  fprintf(' 2.5 percentile \n');
  for i = 1:2
      fprintf('Parameter %1g: %12g \n', [i, S(3,i)])
  end
  fprintf('\n')

  fprintf(' Median \n');
  for i = 1:2
      fprintf('Parameter %1g: %12g \n', [i, S(4,i)])
  end
  fprintf('\n')

  fprintf(' 97.5 percentile \n');
  for i = 1:2
      fprintf('Parameter %1g: %12g \n', [i, S(5,i)])
  end
  fprintf('\n')
Summary information for posterior distribution 
Mean 
Parameter 1:    -0.179049 
Parameter 2:      4.49284 

Standard deviation 
Parameter 1:     0.377056 
Parameter 2:     0.377069 

 2.5 percentile 
Parameter 1:    -0.884899 
Parameter 2:      3.72752 

 Median 
Parameter 1:    -0.191559 
Parameter 2:      4.50619 

 97.5 percentile 
Parameter 1:     0.585931 
Parameter 2:      5.19878 

Generate histogram of samples

Assign number of bins

  nB = 30;
  for j = 1:2

Assign histogram for MCMC samples

      Aa = exp(AA(:, :, j));
      Aa = Aa(:);
      [na,aa] = hist(Aa, nB);
      dela = aa(2) - aa(1);
      nan = na/(sum(na)*dela);
      aaa = (0:dela:aa(end)+5*dela)';

Generate figure

      figure;
      bar(aa, nan)
      hold on
      if j == 1
          plot(aaa,gampdf(aaa,m0a/2,2/m0a))
          hold off
          xlabel('\alpha/1')
          ylabel('p(\alpha)/1')
          title('MCMC samples: \alpha')
          legend('Posterior','Prior')
      elseif j == 2
          plot(aaa,gampdf(aaa,m0d/2,2/(m0d*s0d^2)))
          hold off
          xlabel('\delta/1')
          ylabel('p(\delta)/1')
          title('MCMC samples: \delta')
          legend('Posterior','Prior')
      end
  end

References

[1] 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

MCMCMH Open-Source Software (Release 1.0): Software for measurement uncertainty evaluation

A B Forbes and K. Jagan, Data Science Group, National Physical Laboratory, UK

Copyright NPL Management Limited 2018