Procedure: Recursively update the dynamic emulator mean and variance in the approximation method

Description and Background

This page is concerned with task of emulating a dynamic simulator, as set out in the variant thread on dynamic emulation (ThreadVariantDynamic).

The approximation procedure described in page ProcApproximateIterateSingleStepEmulator recursively defines

\[\mu_{t+1}= \mathrm{E}[ m^*(w_t,a_{t+1},\phi)|f(D),\theta]\]
\[V_{t+1}= \mathrm{E}[ v^*\{(w_t,a_{t+1},\phi),(w_t,a_{t+1},\phi)\}|f(D),\theta] + \mathrm{Var}[m^*(w_t,a_{t+1},\phi)|f(D),\theta]\]

where the expectations and variances are taken with respect to \(w_{t}\), with \(w_{t} \sim N_r(\mu_{t},V_{t})\). Here, we show how to compute \(\mu_{t+1}\) and \(V_{t+1}\) in the case of a single step emulator linear mean and a separable Gaussian covariance function.

To simplify notation, we now omit the constant parameters \(\phi\) from the simulator inputs. (We can think of \(\phi\) as an extra forcing input that is constant over time).

Inputs

  • \(\mu_{t}\) and \(V_{t}\)
  • The single step emulator training inputs \(D\) and outputs \(f(D)\)
  • Emulator covariance function parameters \(B\) and \(\Sigma\) (defined below).

Outputs

  • \(\mu_{t+1}\) and \(V_{t+1}\)

Notation

\(x\): a \(p \times 1\) input vector \(( {x^{(w)}}^T,{x^{(a)}}^T )^T\), with \(x^{(w)}\) the corresponding \(r \times 1\) vector state variable input, and \(x^{(a)}\) the corresponding \((p-r) \times 1\) vector forcing variable input.

\(h(x^{(w)},x^{(a)})= h(x)=(1\, x^T)^T\): the prior mean function.

\(\Sigma c(x,x') =\Sigma \exp\{-(x-x')^TB(x-x')\}\): the prior covariance function, with

\(\Sigma\): the covariance matrix between outputs at the same input \(x\), and

\(B\): a diagonal matrix with \((i,j)^{\mathrm{th}}\) element \(1/\delta_i^2\).

\(B_w\): the upper-left \(r \times r\) submatrix of :math:` B`.

\(B_a\): the lower-right \((p-r) \times (p-r)\) submatrix of \(B\).

\(D\): the set of training data inputs \(x_1,\ldots,n\).

\(c\{(x^{(w)},x^{(a)}),D\}=c(x,D)\): an \(n\times 1\) vector with \(i^{\mathrm{th}}\) element \(c(x,x_i)\).

\(A\): an \(n\times n\) matrix with \((i,j)^{\mathrm{th}}\) element \(c(x_i,x_j)\).

\(H\): an \(n\times (r+1)\) matrix with \(i^{\mathrm{th}}\) row \(h(x_i)^T\).

\(f(D)\): an \(n\times r\) matrix of training outputs with \((i,j)^{\mathrm{th}}\) element the \(j^{\mathrm{th}}\) training output type for the \(i^{\mathrm{th}}\) training input.

\(\hat{\beta}=(H^TA^{-1}H)^{-1}H^TA^{-1}f(D)\).

\(0_{a\times b}\): an \({a\times b}\) matrix of zeros.

Procedure

A series of constants need to be evaluated in the following order. Terms defined at each stage can be evaluated independently of each other, but are expressed in terms of constants defined at preceding stages.

Stage 1: compute \(K_{Vh},\) \(K_{Ec},\) \(K_{Eh},\) \(K_{Ewc}\) and \(K_{Ecc}\)

\[\begin{split}K_{Vh}=\mathrm{Var}[h(w_t,a_{t+1})|f(D),B] = \left(\begin{array}{ccc}0_{1\times 1} & 0_{1\times r} & 0_{1\times (p-r)} \\ 0_{p\times 1} & V_t & 0_{p\times (p-r)} \\ 0_{(p-r)\times 1} & 0_{(p-r)\times r} & 0_{(p-r)\times (p-r)} \end{array}\right),\end{split}\]

\(K_{Ec}=\mathrm{E}[c\{(w_t,a_{t+1}),D\}|f(D),B]\), an \(n \times 1\) vector, with element \(i\) given by

\[\begin{split}|2V_tB_w+I_r|^{-1/2} \exp\{-(a_{t+1} - x_i^{(a)})^TB_a(a_{t+1} - x_i^{(a)})\} \\ \times \exp\{-(\mu_t-x_i^{(w)})^T(2V_t+B_w^{-1})^{-1}(\mu_t - x_i^{(w)})\}\end{split}\]
\[K_{Eh}=\mathrm{E}[h(w_t,a_{t+1})|f(D),B]=(1, \mu_{t}^T,a_{t+1}^T)^T\]

\(K_{Ewc}=\mathrm{E}[w_tc\{(w_t,a_{t+1}),D\}^T|f(D),B]\) , an \(r \times n\) matrix, with column \(i\) given by

\[\begin{split}\mathrm{E}[w_tc(\{w_t,a_{t+1}\},x_i)|f(D),B] = |2V_t B_w + I_r|^{-1/2} \times\exp\{-(a_{t+1}-x_i^{(a)})^TB_a (a_{t+1}-x_i^{(a)}) \} \\ \times \exp\left\{-(\mu_t-x_i^{(w)})^T\left(2V_t+B_W^{-1}\right)^{-1} (\mu_t-x_i^{(w)}) \right\} \times(2B_w+V_t^{-1})^{-1}(2B_wx_i^{(w)} + V_t^{-1}\mu_t).\end{split}\]

\(K_{Ecc}=\mathrm{E}[c\{(w_t,a_{t+1}),D\}c\{(w_t,a_{t+1}),D\}^T |f(D),B]\), an \(n \times n\) matrix, with element \(i,j\) given by

\[\begin{split}\mathrm{E}[c(\{w_t,a_{t+1}\},x_i)c(\{w_t,a_{t+1}\},x_j)|f(D),B] = |4V_t B_w + I_r|^{-1/2} \exp \left\{-\frac{1}{2}(x_i^{(w)}- x_j^{(w)})^TB_w (x_i^{(w)}- x_j^{(w)}) \right\} \\ \times\exp\{-(a_{t+1}-x_i^{(a)})^TB_a (a_{t+1}-x_i^{(a)})- (a_{t+1}-x_j^{(a)})^TB_a (a_{t+1}-x_j^{(a)}) \} \\ \times \exp\left[-\left\{\mu_t-\frac{1}{2}(x_i^{(w)}+ x_j^{(w)}) \right\}^T\left(2V_t+\frac{1}{2}B_W^{-1}\right)^{-1} \left\{\mu_t-\frac{1}{2}(x_i^{(w)}+ x_j^{(w)}) \right\} \right]\end{split}\]

Stage 2: compute \(K_{Cwc},\) \(K_{Ehh},\) and \(K_{Vc}\)

\[K_{Cwc}=\mathrm{Cov}[w_t,c\{(w_t,a_{t+1}),D\}|f(D),B]=K_{Ewc} - \mu_tK_{Ec}^T\]
\[K_{Ehh}=\mathrm{E}[h(w_t,a_{t+1})h(w_t,a_{t+1})^T|f(D),B]=K_{Vh} + K_{Eh}K_{Eh}^T\]
\[K_{Vc}=\mathrm{Var}[c\{(w_t,a_{t+1}),D\}|f(D),B]=K_{Ecc} - K_{Ec}K_{Ec}^T\]

Stage 3: compute \(K_{Chc}\)

\[\begin{split}K_{Chc}=\mathrm{Cov}[h(w_t,a_{t+1}),c\{(w_t,a_{t+1}),D\}|f(D),B] =\left(\begin{array}{c}0_{1\times n} \\ K_{Cwc} \\ 0_{(p-r)\times n} \end{array}\right)\end{split}\]

Stage 4: compute \(K_{Ehc}\) and \(K_{Vm}\)

\[K_{Ehc}=\mathrm{E}[h(w_t,a_{t+1})c\{(w_t,a_{t+1}),D\}^T|f(D),B]=K_{Chc} + K_{Eh}K_{Ec}^T\]
\[\begin{split}K_{Vm}= \mathrm{Var}[m^*(w_t,a_{t+1})|f(D),B] = \hat{\beta}^T K_{Vh}\hat{\beta} +\hat{\beta}^T K_{Chc}A^{-1}(f(D)-H\hat{\beta}) \\ +(f(D)-H\hat{\beta})^TK_{Chc}^T\hat{\beta}+(f(D)- H\hat{\beta})^TA^{-1}K_{Vc}A^{-1}(f(D)-H\hat{\beta})\end{split}\]

Stage 5: compute \(K_{Ev}\)

\[\begin{split}K_{Ev}= \mathrm{E}[v^*\{(w_t,a_{t+1}),(w_t,a_{t+1})\}|f(D),B] = 1 - {\rm tr}[\{A^{-1}-A^{-1}H(H^TA^{-1}H)^{-1}H^TA^{-1}\}K_{Ecc}]\\ + {\rm tr}[(H^TA^{-1}H)^{-1}K_{Ehh} ]-2{\rm tr}[A^{-1}H(H^TA^{-1}H)^{-1}K_{Ehc}]\end{split}\]

Stage 6: compute the procedure outputs \(\mu_{t+1}\) and \(V_{t+1}\)

\[\mu_{t+1} = K_{Eh} \hat{\beta}+K_{Ec}^TA^{-1}(f(D)-H\hat{\beta})\]
\[V_{t+1} = K_{Vm}+K_{Ev}\Sigma\]

Reference

Conti, S., Gosling, J. P., Oakley, J. E. and O’Hagan, A. (2009). Gaussian process emulation of dynamic computer codes. Biometrika 96, 663-676.