Procedure: Sampling the posterior distribution of the correlation lengths

Description and Background

This procedure shows how to draw samples from the posterior distribution of the correlation lengths \(\pi^*_{\delta}(\delta)\) and how to use the drawn samples to make predictions using the emulator. This procedure complements the procedure for building a Gaussian process emulator for the core problem (ProcBuildCoreGP), and represents the formal way of accounting for the uncertainty about the true value of the correlation lengths \(\delta\).

Inputs

  • An emulator as defined in ProcBuildCoreGP, using a linear mean and a weak prior.

Outputs

  • A set of \(s\) samples for the correlation lengths \(\delta\), denoted as \(\tilde{\delta}\).
  • A posterior mean \(\tilde{m}^*(\cdot)\) and covariance function \(\tilde{u}^*(\cdot,\cdot)\), conditioned on the drawn samples \(\tilde{\delta}\).

Procedure

A method for drawing samples from the posterior distribution of \(\delta\) is via the Metropolis-Hastings algorithm. Setting up this algorithm requires an initial estimate of the correlation lengths, and a proposal distribution. An initial estimate can be the value that maximises the posterior distribution, as this is defined in the discussion page for finding the posterior mode of correlation lengths (DiscPostModeDelta). We call this estimate \(\hat{\delta}\).

According to reference [1], a proposal distribution can be

\[\delta^{(i)} \sim {\cal N}(\delta^{(i-1)},V)\]

where \(\delta^{(i-1)}\) is the sample drawn at the \((i-1)\)-th step of the algorithm, and

\[\displaystyle V = -\frac{2.4}{\sqrt{p}}\left( \frac{\partial^2 \pi^*_{\delta}(\hat{\delta})}{\partial \hat{\delta}}\right)^{-1}\]

The Hessian matrix \(\frac{\partial^2 \pi^*_{\delta}(\hat{\delta})}{\partial \hat{\delta}}\) is the same as \(\frac{\partial^2 g(\tau)}{\partial \tau_l \partial \tau_k}\), which was defined in DiscPostModeDelta, after substituting the derivatives \(\partial A / \partial \tau\) with the derivatives \(\partial A / \partial \delta\). The latter are given by

\[\begin{split}\displaystyle \left(\frac{\partial A} {\partial \delta_k}\right)_{i,j} &= A(i,j)\left[\frac{2(x_{k,i}-x_{k,j})^2}{\delta_k^3}\right] \\ \displaystyle \left(\frac{\partial^2 A} {\partial^2 \delta_k}\right)_{i,j} &= A(i,j) \frac{(x_{k,i}-x_{k,j})^2}{\delta_k^4} \left[ \frac{4(x_{k,i}-x_{k,j})^2}{\delta_k^2} - 6 \right]\end{split}\]

and finally

\[\displaystyle \left(\frac{\partial^2 A} {\partial \delta_l\partial \delta_k}\right)_{i,j} = A(i,j) \left[\frac{2(x_{l,i}-x_{l,j})^2}{\delta_l^3}\right] \left[\frac{2(x_{k,i}-x_{k,j})^2}{\delta_k^3}\right]\]

Having defined the initial estimate of \(\delta\) and the proposal distribution \({\cal N}(\delta^{(i-1)},V)\) we can set up the following Metropolis Hastings algorithm

  1. Set \(\delta^{(1)}\) equal to \(\hat{\delta}\)
  2. Add to \(\delta^{(i)}\) a normal variate drawn from \({\cal N}(0,V)\) and call the result \(\delta'\)
  3. Calculate the ratio \(\displaystyle \alpha = \frac{\pi^*_{\delta}(\delta')}{\pi^*_{\delta}(\delta^{(i)})}\)
  4. Draw \(w\) from a uniform distribution in [0,1]
  5. if \(w <\alpha\) set \(\delta^{(i+1)}\) equal to \(\delta'\), else set it equal to \(\delta^{(i)}\)
  6. Repeat steps 2-5 \(s\) times

Finally, if we define as \(m^{*(i)}(x)\) and \(u^{*(i)}(x,x')\) the posterior mean and variance of the emulator using sample \(\delta^{(i)}\), a total estimate of these two quantities, taking into account all the \(s\) samples of \(\delta\) drawn, is given by

\[\displaystyle \tilde{m}^*(x) = \frac{1}{s}\sum_{i=1}^s m^{*(i)}(x)\]

and

\[\displaystyle \tilde{u}^*(x,x') = \frac{1}{s}\sum_{i=1}^s u^{*(i)}(x,x') + \frac{1}{s}\sum_{i=1}^s \left[m^{*(i)}(x) - \tilde{m}^*(x)\right] \left[m^{*(i)}(x') - \tilde{m}^*(x')\right]\]

The procedure for predicting the simulator outputs using more than one hyperparameter sets is described in greater detail in page (ProcPredictGP).

References

  1. Gilks, W.R., Richardson, S. & Spiegelhalter, D.J. (1996). Markov Chain Monte Carlo in Practice. Chapman & Hall.