Procedure: Variance Based Sensitivity Analysis using a GP emulator

Description

This page describes the formulae needed for performing Variance Based Sensitivity Analysis (VBSA). In VBSA we consider the effect on the output \(f(X)\) of a simulator as we vary the inputs \(X\), when the variation of those inputs is described by a (joint) probability distribution \(\omega(X)\). This probability distribution can be interpreted as describing uncertainty about the best or true values for the inputs, or may simply represent the range of input values of interest to us. The principal measures for quantifying the sensitivity of \(f(X)\) to a set of inputs \(X_w\) are the main effect, the sensitivity index and the total effect index. We can also define the interaction between two inputs \(X_i\) and \(X_j\). These are described below.

Mean, main and interaction effects

Let \(w\) denote a subset of the indices from 1 to the number \(p\) of inputs. Let \(X_w\) denote the set of inputs with indices in \(w\). The mean effect of a set of inputs \(X_w\) is the function

\[M_w(x_w) = \textrm{E}[f(X)|x_w].\]

This is a function of \(x_w\) showing how the simulator output for given \(x_w\), when averaged over the uncertain values of all the other inputs, varies with \(x_w\).

The deviation of the mean effect \(M_{\{i\}}(x_{\{i\}})\) of the single input \(X_i\) from the overall mean,

\[I_i(x_w) = \textrm{E}[f(X)|x_i]-\textrm{E}[f(X)].\]

is called the main effect of \(X_i\).

Furthermore, we define the interaction effect between inputs \(X_i\) and \(X_j\) to be

\[I_{\{i,j\}}(x_i,x_j) = \textrm{E}[f(X)|x_i,x_j]-I_i(x_i) - I_j(x_j) - \textrm{E}[f(X)]\]

Sensitivity variance

The sensitivity variance \(V_w\) describes the amount by which the uncertainty in the output is reduced when we are certain about the values of the inputs \(X_w\). That is,

\[V_w = \textrm{Var}[f(X)]-\textrm{E}[\textrm{Var}[f(X)|x_w]]\]

It can also be written in the following equivalent forms

\[V_w = \textrm{Var}[\textrm{E}[f(X)|x_w]] = \textrm{E}[\textrm{E}[f(X)|x_w]^2] - \textrm{E}[f(X)]^2\]

Total effect variance

The total effect variance \(V_{Tw}\) is the expected amount of uncertainty in the model output that would be left if we removed the uncertainty in all the inputs except for the inputs with indices \(w\). The total effect variance is given by

\[V_{Tw} = \textrm{E}[\textrm{Var}[f(X)|x_{\bar{w}}]]\]

where \(\bar{w}\) denotes the set of all indices not in \(w\), and hence \(x_{\bar{w}}\) means all the inputs except for those in \(x_w\).

\(V_{Tw}\) can also be written as

\[V_{Tw} = \textrm{Var}[f(X)] - V_{\bar{w}}.\]

In order to define the above quantities, it is necessary to specify a full probability distribution for the inputs. This clearly demands a good understanding of probability and its use to express personal degrees of belief. However, this specification of uncertainty often also requires interaction with relevant experts whose knowledge is being used to specify values for individual inputs. There is a considerable literature on the elicitation of expert knowledge and uncertainty in probabilistic form, and some references are given at the end of the procedure page for uncertainty analysis (ProcUAGP) page.

We assume here that a Gaussian process (GP) emulator has been built according to the procedure page ProcBuildCoreGP, and that we are only emulating a single output. Note that ProcBuildCoreGP gives procedures for deriving emulators in a number of different forms, and we consider here only the “linear mean and weak prior” case where the GP has a linear mean function, weak prior information is specified on the hyperparameters \(\beta\) and \(\sigma^2\) and the emulator is derived with a single point estimate for the hyperparameters \(\delta\).

The procedure here computes the expected values (with respect to the emulator) of the above quantities.

Inputs

• An emulator as defined in ProcBuildCoreGP
• A distribution \(\omega(\cdot)\) for the uncertain inputs
• A set \(w\) of indices for inputs whose average effect or sensitivity indices are to be computed, or a pair {i,j} of indices defining an interaction effect to be computed
• Values \(x_w\) for the inputs \(X_w\), or similarly for \(X_i,X_j\)

Outputs

• \(\textrm{E}^*[M_w(x_w)]\)
• \(\textrm{E}^*[I_{\{i,j\}}(x_i,x_j)]\)
• \(\textrm{E}^*[V_w]\)
• \(\textrm{E}^*[V_{Tw}]\)

where \(\textrm{E}^*[\cdot]\) denotes an expectation taken with respect to the emulator uncertainty, i.e. a posterior mean.

Procedure

In this section we provide the formulae for the calculation of the posterior means of \(M_w(x_w)\), \(I_{\{i,j\}}(x_i,x_j)\), \(V_w\) and \(V_{Tw}\). These are given as a function of a number of integral forms, which are denoted as \(U_w,P_w,S_w,Q_w,R_w\) and \(T_w\). The exact expressions for these forms depend on the distribution of the inputs \(\omega(\cdot)\), the correlation function \(c(.,.)\) and the regression function \(h(\cdot)\). In the following section, we give expressions for the above integral forms for the general and two special cases.

Calculation of \(\textrm{E}^*[M_w(x_w)]\)

\[\textrm{E}^*[M_w(x_w)] = R_w\hat{\beta} + T_w e,\]

where \(e = A^{-1}(f(D)-H\hat{\beta})\) and \(\hat{\beta}, A, f(D)\) and \(H\) are defined in ProcBuildCoreGP.

For the main effect of \(X_i\) the posterior mean is

\[\textrm{E}^*[I_i(x_i)] = \{R_{\{i\}} - R\}\hat{\beta} + \{T_{\{i\}}-T\}e.\]

It is important to note here that both \(R_w\) and \(T_w\) are functions of \(x_w\). The dependence on \(x_w\) has been suppressed here for notational simplicity.

Calculation of \(\textrm{E}^*[I_{\{i,j\}}(x_i,x_j)]\)

\[\begin{split}\begin{array}{r l} \textrm{E}^*[I_{\{i,j\}}(x_i,x_j)] = & \{R_{\{i,j\}} - R_{\{i\}} - R_{\{j\}} - R\}\hat{\beta} \\ + & \{T_{\{i,j\}} - T_{\{i\}} - T_{\{j\}} - T\}e \end{array}\end{split}\]

where \(R_{\{i,j\}}\) and \(R_{\{i\}}\), for instance, are special cases of \(R_w\) when the set \(w\) of indices comprises the two elements \(w=\{i,j\}\) or the single element \(w=\{i\}\). Remember also that these will be functions of \(x_{\{i,j\}}=(x_i,x_j)\) and \(x_{\{i\}}=x_i\) respectively.

Calculation of \(\textrm{E}^*[V_w]\)

We write the posterior mean of \(V_w\) as

\[\textrm{E}^*[V_w] = \textrm{E}^*[\textrm{E}[\textrm{E}[f(X)|x_w]^2]] - \textrm{E}^*[\textrm{E}[f(X)]^2]\]

The first term is

\[\begin{split}\begin{array}{r l} \textrm{E}^*[\textrm{E}[\textrm{E}[f(X)|x_w]^2]] & = \hat{\sigma}^2[U_w-\mathrm{tr}(A^{-1}P_w) + \mathrm{tr}\{W(Q_w-S_w A^{-1} H \\ & \qquad \qquad - H^{\mathrm{T}}A^{-1}S_w^{\mathrm{T}} + H^{\mathrm{T}}A^{-1}P_w A^{-1}H)\}] \\ & \quad + e^{\mathrm{T}}P_we + 2\hat{\beta}^{\mathrm{T}}S_we + \hat{\beta}^{\mathrm{T}}Q_w\hat{\beta} \end{array}\end{split}\]

where \(\hat\sigma^2\) is defined in ProcBuildCoreGP.

The second term is

\[\begin{split}\begin{array}{r l} \textrm{E}^*[\textrm{E}[f(X)]^2] & = \hat{\sigma}^2[U-TA^{-1}T^{\mathrm{T}} +\{R - TA^{-1}H\}W\{R - TA^{-1}H\}^\mathrm{T}] \\ & \quad + \left(R\hat{\beta}+Te\right)^2 \end{array}\end{split}\]

with \(W = (H^{\mathrm{T}}A^{-1}H)^{-1}\).

Calculation of \(\textrm{E}^*[V_{Tw}]\)

\(\textrm{E}^*[V_{Tw}]\) can be calculated via the sensitivity variance \(V_{\bar{w}}\) using the relation

\[\textrm{E}^*[V_{Tw}] = \textrm{E}^*[\textrm{Var}[f(X)]] - \textrm{E}^*[V_{\bar{w}}]\]

with

\[\textrm{E}^*[\textrm{Var}[f(X)]] = \hat{\sigma}^2[U-TA^{-1}T^{\mathrm{T}} +\{R - TA^{-1}H\}W\{R - TA^{-1}H\}^\mathrm{T}]\]

Dimensions

Before presenting the integral forms that appear in the above expressions, we give the dimensions of all the involved quantities in the table below. We assume that we have \(n\) observations, \(p\) inputs and \(q\) regression functions. The terms in the left column are either described in ProcBuildCoreGP or they are shorthands \((e, W)\). The terms in the right hand side column are the integral forms, which will be presented in the following section.

Symbol	Dimension	Symbol	Dimension
\(\hat{\sigma}^2\)	\(1 \times 1\)	\(U_w\)	\(1 \times 1\)
\(\hat{\beta}\)	\(q \times 1\)	\(P_w\)	\(n \times n\)
\(f\)	\(n \times 1\)	\(S_w\)	\(q \times n\)
\(H\)	\(n \times q\)	\(Q_w\)	\(q \times q\)
\(A\)	\(n \times n\)	\(R_w\)	\(1 \times q\)
\(e\)	\(n \times 1\)	\(T_w\)	\(1 \times n\)
\(W\)	\(q \times q\)

The integral forms

General case

When no assumptions are made about the distribution of the inputs, the correlation and the regression functions we have general expressions for the \(U_w, P_w, S_w, Q_w, R_w, T_w\) terms. These are

\[U_w = \int_{{\cal X}_w}\int_{{\cal X}_{\bar{w}}}\int_{{\cal X}_{\bar{w}}} c(x,x^*)\omega(x_{\bar{w}}|x_w)\omega(x'_{\bar{w}}|x_w)\omega(x_w) \mathrm{d} x_{\bar{w}} \mathrm{d} x'_{\bar{w}} \mathrm{d} x_w\]

\[P_w = \int_{{\cal X}_w}\int_{{\cal X}_{\bar{w}}}\int_{{\cal X}_{\bar{w}}} t(x)t(x^*)^{\mathrm{T}} \omega(x_{\bar{w}}|x_w)\omega(x'_{\bar{w}}|x_w)\omega(x_w) \mathrm{d} x_{\bar{w}} \mathrm{d} x'_{\bar{w}} \mathrm{d} x_w\]

\[S_w = \int_{{\cal X}_w}\int_{{\cal X}_{\bar{w}}}\int_{{\cal X}_{\bar{w}}} h(x)t(x^*)^{\mathrm{T}} \omega(x_{\bar{w}}|x_w)\omega(x'_{\bar{w}}|x_w)\omega(x_w) \mathrm{d} x_{\bar{w}} \mathrm{d} x'_{\bar{w}} \mathrm{d} x_w\]

\[Q_w = \int_{{\cal X}_w}\int_{{\cal X}_{\bar{w}}}\int_{{\cal X}_{\bar{w}}} h(x)h(x^*)^{\mathrm{T}} \omega(x_{\bar{w}}|x_w)\omega(x'_{\bar{w}}|x_w)\omega(x_w) \mathrm{d} x_{\bar{w}} \mathrm{d} x'_{\bar{w}} \mathrm{d} x_w\]

\[R_w = \int_{{\cal X}_{\bar{w}}} h(x)^{\mathrm{T}}\omega(x_{\bar{w}}|x_w) \mathrm{d} x_{\bar{w}}\]

\[T_w = \int_{{\cal X}_{\bar{w}}} t(x)^{\mathrm{T}}\omega(x_{\bar{w}}|x_w) \mathrm{d} x_{\bar{w}}\]

Here, \(x_{\bar{w}}\) and \(x'_{\bar{w}}\) denote two different realisations of \(x_{\bar{w}}\). \(x^*\) is a vector with elements made up of \(x'_{\bar{w}}\) and \(x_w\) in the same way as \(x\) is composed of \(x_{\bar{w}}\) and \(x_w\). Remember also that \(R_w\) and \(T_w\) are functions of \(x_w\).

\(h(x)\) is described in the alternatives page on emulator prior mean function (AltMeanFunction). \(c(.,.)\) is the correlation function discussed in the alternatives page on emulator prior correlation function (AltCorrelationFunction). Also \(t(x) = c(D,x)\), as introduced in ProcBuildCoreGP.

\(\omega(x_w)\) is the joint distribution of the \(x_w\) inputs and \(\omega(x_{\bar{w}}|x_w)\) is the conditional distribution of \(x_{\bar{w}}\) when the values of \(x_{w}\) are known.

Finally, when one of the above integral forms appears without a subscript (e.g. \(U\)), it is implied that the set \(w\) is empty.

Special case 1

We now show derive closed form expressions for the above integrals when we make specific choices about \(\omega(\cdot)\) \(c(\cdot,\cdot)\) and \(h(\cdot)\). We first assume that \(\omega(\cdot)\) is a normal distribution given by

\[\omega(x) = \frac{1}{(2\pi)^{d/2}|B|^{-1/2}}\exp\left[-\frac{1}{2}(x-m)^T B (x-m)\right]\]

We also assume the generalised Gaussian correlation function with nugget (see AltCorrelationFunction)

\[c(x,x') = \nu I_{x=x'} + (1-\nu)\exp\{-(x-x')^T C (x-x')\}\]

where \(I_{x=x'}\) equals 1 if \(x=x'\) but is otherwise zero, and \(\nu\) represents the nugget term. The case of a generalised Gaussian correlation function without nugget is simply obtained by setting \(\nu=0\).

We let both \(B,C\) be general positive definite matrices, partitioned as

\[\begin{split}B = \left[ \begin{array}{cc} B_{ww} & B_{w\bar{w}} \\ B_{\bar{w}w} & B_{\bar{w}\bar{w}} \end{array} \right], \\ C = \left[ \begin{array}{cc} C_{ww} & C_{w\bar{w}} \\ C_{\bar{w}w} & C_{\bar{w}\bar{w}} \end{array} \right]\end{split}\]

Finally, we do not make any particular assumption for \(h(x)\).

We now give the expressions for each of the integrals

---

\(U_w\) is the scalar

\[U_w = (1-\nu)\frac{|B|^{1/2}|B_{\bar{w}\bar{w}}|^{1/2}}{|F|^{1/2}}\]

with

\[\begin{split}F = \left[ \begin{array}{ccc} B_{ww} + B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w} & B_{w\bar{w}} & B_{w\bar{w}} \\ B_{\bar{w}w} & 2C_{\bar{w}\bar{w}} + B_{\bar{w}\bar{w}} & -2C_{\bar{w}\bar{w}} \\ B_{\bar{w}w} & -2C_{\bar{w}\bar{w}} & 2C_{\bar{w}\bar{w}} + B_{\bar{w}\bar{w}} \end{array} \right]\end{split}\]

\(U\) is the special case when \(w\) is the empty set. The exact formula for \(U\) is given in ProcUAGP.

---

\(P_w\) is an \(n \times n\) matrix, whose \((k,l)^{\mathrm{th}}\) entry is

\[P_w(k,l) = (1-\nu)^2\frac{|B|^{1/2}|B_{\bar{w}\bar{w}}|^{1/2}}{|F|^{1/2}} \exp\left\{-\frac{1}{2}\left[ r - g^{\mathrm{T}}F^{-1}g \right]\right\}\]

with

\[\begin{split}F = \left[ \begin{array}{ccc} 4C_{ww} + B_{ww} + B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w} & 2C_{w\bar{w}} + B_{w\bar{w}}& 2C_{w\bar{w}} + B_{w\bar{w}} \\ 2C_{\bar{w}w} + B_{\bar{w}w} & 2C_{\bar{w}\bar{w}} + B_{\bar{w}\bar{w}}& 0 \\ 2C_{\bar{w}w} + B_{\bar{w}w} & 0 & 2C_{\bar{w}\bar{w}} + B_{\bar{w}\bar{w}} \end{array} \right]\end{split}\]

and

\[\begin{split}g = \left[ \begin{array}{l} 2C_{ww}(x_{w,k}+x_{w,l} - 2m_w) + 2C_{w\bar{w}}(x_{\bar{w},k}+x_{\bar{w},l} - 2m_{\bar{w}}) \\ 2C_{\bar{w}w}(x_{w,k} - m_w) + 2C_{\bar{w}\bar{w}}(x_{\bar{w},k} - m_{\bar{w}}) \\ 2C_{\bar{w}w}(x_{w,l} - m_w) + 2C_{\bar{w}\bar{w}}(x_{\bar{w},l} - m_{\bar{w}}) \end{array} \right]\end{split}\]

and

\[r = (x_k - m)^{\mathrm{T}}2C(x_k - m) + (x_l - m)^{\mathrm{T}}2C(x_l - m)\]

\(P\) is a special case of \(P_w\) when \(w\) is the empty set, and reduces to

\[P=T^T T\]

---

\(S_w\) is an \(q \times n\) matrix, whose \((k,l)^{\mathrm{th}}\) entry is

\[S_w(k,l) = (1-\nu)\frac{|B|^{1/2}|B_{\bar{w}\bar{w}}|^{1/2}}{|F|^{1/2}} \exp \left\{-\frac{1}{2}\left[r - g^{\mathrm{T}}F^{-1}g\right]\right\} \textrm{E}_*[h_k(x)]\]

with

\[\begin{split}F = \left[ \begin{array}{ccc} 2C_{ww}+B_{ww} + B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w} & B_{w\bar{w}} & 2C_{w\bar{w}}+B_{w\bar{w}} \\ B_{\bar{w}w} & B_{\bar{w}\bar{w}} & 0 \\ 2C_{\bar{w}w}+B_{\bar{w}w} & 0 & 2C_{\bar{w}\bar{w}}+B_{\bar{w}\bar{w}} \end{array} \right]\end{split}\]

and

\[\begin{split}g = \left[\begin{array}{ccc} 2C_{ww} & 0 & 2C_{w\bar{w}} \\ 0 & 0 & 0 \\ 2C_{\bar{w}w} & 0 & 2C_{\bar{w}\bar{w}} \end{array}\right] \left[ \begin{array}{c} x_{w,l} - m_w \\ 0 \\ x_{\bar{w},l} - m_{\bar{w}} \end{array} \right]\end{split}\]

and

\[\begin{split}r = \left[ \begin{array}{c} x_{w,l} - m_w \\ 0 \\ x_{\bar{w},l} - m_{\bar{w}} \end{array} \right]^{\mathrm{T}} \left[\begin{array}{ccc} 2C_{ww} & 0 & 2C_{w\bar{w}} \\ 0 & 0 & 0 \\ 2C_{\bar{w}w} & 0 & 2C_{\bar{w}\bar{w}} \end{array}\right] \left[ \begin{array}{c} x_{w,l} - m_w \\ 0 \\ x_{\bar{w},l} - m_{\bar{w}} \end{array} \right]\end{split}\]

The expectation \(\textrm{E}_z[\cdot]\) is w.r.t. the normal distribution \({\cal{N}}(m + F^{-1}g,F^{-1})\). Also \(h_k(x)\) is the \(k\)-th element of \(h(x)\).

\(S\) is a special case of \(S_w\) when \(w\) is the empty set, and reduces to

\(S=R^T T\)

---

\(Q_w\) is a \(q \times q\) matrix, whose \((k,l)^{\mathrm{th}}\) entry is

\[Q_w(k,l) = \frac{|B|^{1/2}|B_{\bar{w}\bar{w}}|^{1/2}}{|F|^{1/2}} \textrm{E}_*[h_k(x)h_l^{\mathrm{T}}(x)]\]

where the expectation \(\textrm{E}_*[\cdot]\) is w.r.t. the normal distribution \({\cal{N}}([m_w,m_{\bar{w}}]^{\mathrm{T}},F^{-1})\)

with

\[\begin{split}F = \left[ \begin{array}{cc} B_{ww} + B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w} & B_{w\bar{w}} \\ B_{\bar{w}w} & B_{\bar{w}\bar{w}} \end{array} \right]\end{split}\]

\(Q\) is a special case of \(Q_w\) when \(w\) is the empty set, and reduces to

\[Q=R^T R\]

---

\(R_w\) is the \(1 \times q\) vector with elements the mean of the elements of \(h(x)\), w.r.t. \(\omega(x_{\bar{w}}|x_w)\), i.e.,

\[R_w = \int_{{\cal X}_{\bar{w}}} h(x)^{\mathrm{T}}\omega(x_{\bar{w}}|x_w) \mathrm{d} x_{\bar{w}}\]

and is a function of \(x_w\). \(R\) is a special case of \(R_w\), when \(w\) is the empty set. The formula for \(R\) is given in ProcUAGP.

---

\(T_w\) is an \(1 \times n\) vector, whose \(k^{\mathrm{th}}\) entry is

\[T_w(k) = (1-\nu) \frac{|B_{\bar{w}\bar{w}}|^{1/2}}{|2C_{\bar{w}\bar{w}}+B_{\bar{w}\bar{w}}|^{1/2}} \exp\left\{-\frac{1}{2}\left[F^{'} +r-g^TF^{-1}g\right]\right\}\]

with

\[\begin{split}F^{'} &= (x_w-m_w - (F^{-1}g)_w)^T \\ & \big[2C_{ww}+B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w} - (2C_{w\bar{w}} +B_{w\bar{w}})(2C_{\bar{w}\bar{w}} +B_{\bar{w}\bar{w}})^{-1}(2C_{\bar{w}w} +B_{\bar{w}w})\big] \\ & (x_w-m_w - (F^{-1}g)_w) \\ F &= \left[ \begin{array}{cc} 2C_{ww} + B_{w\bar{w}}B_{\bar{w}\bar{w}}^{-1}B_{\bar{w}w}& 2C_{w\bar{w}} + B_{w\bar{w}} \\ 2C_{\bar{w}w} + B_{\bar{w}w}&2C_{\bar{w}\bar{w}} + B_{\bar{w}\bar{w}} \end{array} \right] \\ g &= 2C(x_k-m) \\ r &= (x_k-m)^T 2C(x_k-m)\end{split}\]

\((F^{-1}g)_w\) is the part of the \(F^{-1}g\) vector that corresponds to the indices \(w\). According to the above formulation, these are the first \(\#(w)\) indices, where \(\#(w)\) is the number of indices contained in \(w\).

\(T\) is a special case of \(T_w\), when \(w\) is the empty set. The formula for \(T\) is given in ProcUAGP.

Special case 2

In this special case, we further assume that the matrices \(B,C\) are diagonal. We also consider a special form for the vector \(h(x)\), which is the linear function described in AltMeanFunction

\[h(x) = [1,x^{\mathrm{T}}]^{\mathrm{T}}\]

We now present the form of the integrals under the new assumptions.

---

\(U_w\) is the scalar

\[U_w = (1-\nu)\prod_{i\in \bar{w}} \left(\frac{B_{ii}}{B_{ii} + 2(2C_{ii})}\right)^{1/2}\]

Again, \(U\) is the special case when \(w\) is the empty set, and its exact formula is given in ProcUAGP.

---

\(P_w\) is an \(n \times n\) matrix, whose \((k,l)^{\mathrm{th}}\) entry is

\[\begin{split}\begin{array}{r l} P_w(k,l) = &(1-\nu)^2\prod_{i\in {\bar{w}}} \frac{B_{ii}}{2C_{ii}+B_{ii}} \exp\left\{-\frac{1}{2}\frac{2C_{ii} B_{ii}}{2C_{ii}+B_{ii}} \left[(x_{i,k}-m_i)^2 + (x_{i,l}-m_i)^2\right]\right\} \\ & \prod_{i\in w} \left(\frac{B_{ii}}{4C_{ii}+B_{ii}}\right)^{1/2} \exp\left\{-\frac{1}{2}\frac{1}{4C_{ii}+B_{ii}}\right. \\ & \left[4C_{ii}^2(x_{i,k}-x_{i,l})^2 + 2C_{ii} B_{ii}\left\{(x_{i,k}-m_i)^2 + (x_{i,l}-m_i)^2\right\}\right]\Big\} \end{array}\end{split}\]

where the double indexed \(x_{i,k}\) denotes the \(i^{\mathrm{th}}\) input of the \(k^{\mathrm{th}}\) training data.

\(P\) is a special case of \(P_w\) when \(w\) is the empty set, and reduces to

\(P=T^T T\)

---

\(S_w\) is an \(q \times n\) matrix whose \((k,l)^{\mathrm{th}}\) entry is

\[\begin{split}\begin{array}{rl} S_w{(k,l)} = &(1-\nu) \textrm{E}_*[h_k(x)] \\ &\quad\prod_{i\in \{w,\bar{w}\}}\frac{B_{ii}^{1/2}} {(2C_{ii} + B_{ii})^{1/2}} \exp \left\{-\frac{1}{2}\left[\frac{2C_{ii}B_{ii}}{2C_{ii}+B_{ii}}(x_{i,l} - m_i)^2\right]\right\} \end{array}\end{split}\]

For the expectation we have

\[\begin{split}\begin{array}{ll} \textrm{E}_*[h_1(x)] = 1 & \\ \textrm{E}_*[h_{j+1}(x)] = m_j & \qquad \mathrm{if} \quad j \in \bar{w} \\ \textrm{E}_*[h_{j+1}(x)] = \frac{2C_{jj}x_{j,l} + B_{jj}m_j}{2C_{jj} + B_{jj}}& \qquad \mathrm{if} \quad j \in w \end{array}\end{split}\]

\(S\) is a special case of \(S_w\) when \(w\) is the empty set, and reduces to

\(S=R^T T\)

---

\(Q_w\) is a \(q \times q\) matrix. If we assume that its \(q\) indices have the labels \([1,\bar{w},w]\), then,

\[\begin{split}Q_w = \left [ \begin{array}{ccc} 1 & m_{\bar w}^T & m_w^T \\ m_{\bar w} & m_{\bar w}m_{\bar w}^T & m_{\bar w}m_w^T \\ m_w & m_wm_{\bar w}^T & m_wm_w^T + B_{ww}^{-1} \end{array} \right ]\end{split}\]

\(Q\) is a special case of \(Q_w\), when \(w\) is the empty set, and reduces to

\(Q=R^T R\)

---

\(R_w\) is a \(1 \times q\) vector. If we assume that its \(q\) indices have the labels \([1,\bar{w},w]\), then,

\[R_w = [1,m_{\bar{w}}^T,x_w^T]\]

\(R\) is a special case of \(R_w\), when \(w\) is the empty set. The formula for \(R\) is given in ProcUAGP.

---

\(T_w\) is an \(1 \times n\) vector, whose \(k^{\mathrm{th}}\) entry is

\[\begin{split}\begin{array}{r l} \displaystyle T_w(k) = (1-\nu) \prod_{i\in \{\bar{w}\}} \frac{\displaystyle B_{ii}^{1/2}}{\displaystyle (2C_{ii}+B_{ii})^{1/2}} & \displaystyle \exp\left\{-\frac{1}{2}\frac{\displaystyle 2C_{ii}B_{ii}}{\displaystyle 2C_{ii}+B_{ii}} \left(x_{i,k}-m_i\right)^2\right\} \\ & \displaystyle \exp\left\{-\frac{1}{2}(x_w-x_{w,k})^T 2C_{ww} (x_w-x_{w,k})\right\} \end{array}\end{split}\]

Recall that \(x_w\) denotes the fixed values for the inputs \(X_w\), upon which the measures \(M_w\), \(V_w\) and \(V_{Tw}\) are conditioned. On the other hand, \(x_{w,k}\) represents the \(w\) inputs of the \(k\)th design points.

\(T\) is a special case of \(T_w\), when \(w\) is the empty set. The formula for \(T\) is given in ProcUAGP.

References

The principal reference for these procedures is

Oakley, J.E., O’Hagan, A., (2004), Probabilistic Sensitivity Analysis of Complex Models: a Bayesian Approach, J.R. Statist. Soc. B, 66, Part 3, pp.751-769.

The above paper does not explicitly consider the case of a non-zero nugget. The calculations of \({\rm E}^*[V_w]\) and \({\rm E}^*[V_{Tw}]\) produce results that are scaled by \((1-\nu)\), and in general \((1- \nu)\sigma^2\) is the maximum variance reduction achievable because the nugget \(\nu\) represents noise that we cannot learn about by reducing uncertainty about \(X\). See the discussion page on the nugget effects in sensitivity analysis (DiscUANugget) for more details on this point.

Ongoing work

We intend to provide procedures relaxing the assumption of the “linear mean and weak prior” case of ProcBuildCoreGP as part of the ongoing development of the toolkit. We also intend to provide procedures for computing posterior variances of the various measures.