Emulation Development with ExeUQ Package

Build the emulator object

The emulator requires a mean function, a form for the covariance function, a method, and priors for the parameters. The EMULATE.gpstan function construct emulator in a fully Bayesian way and it contains the posterior samples of the parameters. The meanResponse is our linear function computed previously, additionalVariables asks which inputs, in addition to those fitted in the meanResponse, are more active, i.e. have more effect on the model response. FastVersion builds all of the required covariance matrices/cholesky factors at the MAP (Maximum a posteriori estimation), so that fast predictions and diagnostics are possible.

myem.gp.cnrm1 = EMULATE.gpstan(meanResponse=myem.lm.cnrm, tData=tData, additionalVariables=cands[-10], FastVersion = TRUE)

Diagnostics for emulator object

We could easily perform Leave One Put diagnostic plots against each of the model inputs. The predictions and two standard deviation prediction intervals are in black. The model values are in green if they lie within two standard deviation prediction intrevals, or red otherwise.

tLOOs <- LOO.plot(StanEmulator = myem.gp.cnrm1, ParamNames=cands[-10])

To predict with an emulator we use the function EMULATOR.gpstan

validation <- ValidationStan(tData.valid, myem.gp.cnrm1)

## 
## SAMPLING FOR MODEL 'PredictGen' NOW (CHAIN 1).
## Iteration: 1 / 1 [100%]  (Sampling)
## 
##  Elapsed Time: 0 seconds (Warm-up)
##                3.50063 seconds (Sampling)
##                3.50063 seconds (Total)

To identify the problems with the emulator performance we could perform the traditional standardized Leave One Out (LOO) diagnostics fitted to our design data. We need to load key functions first.

source('BuildEmulator/BuildEmulatorNSt.R')

## Loading required package: tgp

## Loading required package: CGP

We use the CalStError, where LOO is the output from LOO.plot function and true.y is the model values at the training (Design) set. From the standardized errors plots we observe the heteroscedasticity in response to ALFX input, this clearly indicates the failure of our stationary (standard) GP emulator.

std.err <- CalStError(LOOs = tLOOs, true.y = tData$SCMdata)
par(mfrow = c(3, 3), mar = c(4, 4, 2, 1))
for(i in 1:(length(cands)-1)) plot(tData[, cands[i]], std.err, pch = 20, xlab = cands[i], ylab = 'std. error')

We are interested in capturing the distinct regional behaviour by observing the standardized errors, i.e. the standardized errors are close to zero for ALFX less than zero and rapidly increase for increasing values of ALFX.

We develop a mixture model that clusters the input points according to the values of standardized errors. For each input point under consideration in X, we derive a vector of probabilities indicating the dominant local behaviour around the point. For this purpose we use function MIXTURE.design, where y.std is the vector of standardized errors and L is the number of input regions of distinct model behaviour.

myem.mixture <- MIXTURE.design(y.std = std.err, X = tData[, cands[-10]], L=2)
mixture.cnrm <- data.frame(cbind(std.err, myem.mixture$MixtureMat, tData[, cands[-10]]))
names(mixture.cnrm) <- c('std.err', 'red', 'blue', cands[-10])
head(mixture.cnrm)

##      std.err         red       blue        ALFX      TENTRX      TENTR
## 1  0.4727585 0.319754099 0.68024590  0.30001137 -0.97069174 -0.1928708
## 2  0.8383965 0.945288062 0.05471194 -0.89968906  0.01740040 -0.4774003
## 3  0.2285943 0.053898109 0.94610189  0.74270569 -0.22388121  0.2830404
## 4 -1.3469797 0.006053914 0.99394609  0.89906039  0.80108132  0.3328454
## 5 -0.1416202 0.427629429 0.57237057 -0.02982919 -0.03334369 -0.6465328
## 6  0.5703032 0.137078422 0.86292158  0.47187173  0.51990493 -0.8869878
##         RKDX        RKDN         VVX        VVN     GAMAP1   REFLCAPE
## 1  0.1205414 -0.64224285 -0.12140496  0.1728073 -0.8185718 -0.9155161
## 2  0.8861751 -0.69197355  0.09457564  0.7393054 -0.7416000  0.6464448
## 3 -0.9669966  0.48148095 -0.09919265 -0.9016202  0.4738241  0.3092267
## 4  0.6318558  0.26913381 -0.62617477 -0.1692253  0.7816716  0.7815247
## 5 -0.6623497  0.56450136  0.84361558 -0.6328080 -0.4242097  0.8602944
## 6 -0.3549518 -0.08295647 -0.16335517  0.6934439 -0.5941942  0.7033497

ggplot(mixture.cnrm, aes(x = ALFX, y = std.err, col = rgb(red = red, blue = blue, green = 0), ymin = -3.5, ymax = 3.5)) + geom_point(size = 2) + scale_color_identity()

We produce the plot of coloured standardized errors against ALFX, where the deep blue colour corresponds to the higher probability of a point being allocated to region 1 (region of high error variability) and the deep red colour corresponds to the higher probability of a point being allocated to region 2 (region of low error variability).

We partitioned the input space into L input regions and for each input region we defined stationary covariance function. For our standard GP emulator we have the following covariance function: \[k(x, x';\sigma^{2}, \delta)= \sigma^{2}r(x, x'; \delta) + \tau^{2}\mathbb{1}\{x=x'\}\]

For our nonstationary GP emulator we have the modified covariance function: \[ k(x, x'; \boldsymbol{\sigma}^2, \boldsymbol{\delta}) = \sum_{l=1}^L \lambda_l(x)\lambda_l(x')\sigma_l^2r_l(x, x';\delta_l) + \mathbb{1}\{x=x'\}\sum_{l=1}^Lz_l(x)z_l(x')\tau_l^2\] where \(z_l(x) = \mathbb{1}(\lambda_l(x)=\max_k(\lambda_k(x)))\).

Such specification allows us to have non-zero covariance between the points from different regions. The function myem.mixture produces matrix myem.mixture$MixtureMat, that contains \(\lambda_l(x), l=1, \dots, L\) for the design matrix.

Build the nonstationry emulator object

The EMULATE.gpstanNSt contains exactly the same inputs as the EMULATE.gpstan, except it requires mixtureComp object found previously.

myem.gp.cnrm.nst1 <- EMULATE.gpstanNSt(meanResponse=myem.lm.cnrm, tData=tData, additionalVariables=cands[-10],
                                       FastVersion = TRUE, mixtureComp = myem.mixture)

Diagnostics for nonstationary emulator object

We could easily perform Leave One Put diagnostic plots against each of the model inputs for nonstationary GP emulator. The predictions and two standard deviation prediction intervals are in black. The model values are in green if they lie within two standard deviation prediction intrevals, or red otherwise.

tLOOs.nst <- LOO.plot.NSt(StanEmulator = myem.gp.cnrm.nst1, ParamNames=cands[-10])

To predict with our nonstationary emulator we are required to compute the mixture components \(\lambda_l(x), l=1, \dots, L\) for our validation points using function MIXTURE.predict and then use the function ValidationStanNSt:

mixture.predict1 <- MIXTURE.predict(x = tData.valid[, cands[-10]], mixtureComp = myem.mixture)

## 
## SAMPLING FOR MODEL 'MixtureModPredict' NOW (CHAIN 1).
## Iteration: 1 / 1 [100%]  (Sampling)
## 
##  Elapsed Time: 0 seconds (Warm-up)
##                0.05674 seconds (Sampling)
##                0.05674 seconds (Total)

validation.nst <- ValidationStanNSt(tData.valid, Emulator=myem.gp.cnrm.nst1, 
                                    mixturePredict=mixture.predict1)

## 
## SAMPLING FOR MODEL 'PredictGenNGP' NOW (CHAIN 1).
## Iteration: 1 / 1 [100%]  (Sampling)
## 
##  Elapsed Time: 0 seconds (Warm-up)
##                45.3595 seconds (Sampling)
##                45.3595 seconds (Total)

Let’s compare the performance of stationary GP emulator against nonstationary GP emulator:

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
ValidPlot(validation, tData.valid[, cands[-10]], tData.valid$SCMdata, axis = 1, 
            heading = "st-GP", xrange = c(-1, 1), interval = range(validation, validation.nst), 
            ParamNames = names(tData.valid))
ValidPlot(validation.nst, tData.valid[, cands[-10]], tData.valid$SCMdata, axis = 1, 
            heading = "nst-GP", xrange = c(-1, 1), interval = range(validation, validation.nst), 
            ParamNames = names(tData.valid))