Example: A one dimensional emulator¶
Model (simulator) description¶
In this page we present an example of fitting an emulator to a \(p=1\) dimensional simulator (model). The simulator we use is surfebm. Surfebm is an energy balance model of the Earth’s climate. The state variables are upper ocean temperatures averaged around the globe at constant latitudes. There are 18 boxes from the North to South poles. In addition there is a single deep ocean box. If the upper ocean boxes have a temperature below 0 it freezes. Water and ice have different albedos and there is an overturning circulation moving heat from the Arctic to the Antarctic. In this example we will keep all inputs fixed, apart from one, the solar constant. For simplicity in this example we only look at a single output, the mean upper ocean temperature. Figure 1 shows the mean temperature over the input range of the solar constant.
Figure 1: The simulator’s output and the 6 design points
Design¶
Design is the selection of the input points at which the simulator is to be run. There are several design options that can be used, as described in the alternatives page on training sample design for the core problem (AltCoreDesign). For this example, we will use an Optimised Latin Hypercube Design, which for one dimension is a set of equidistant points on the space of the input variable.
We select \(n=6\) design points, which are shown in Figure 1 as circles. These points are
We scale our design points so that they lie in \([0,1]\). That is,
The output of the model at these points is
Gaussian Process setup¶
In setting up the Gaussian process, we need to define the mean and the covariance function. For the mean function, we choose the linear form described in the alternatives page for emulator prior mean function (AltMeanFunction), which is \(h(x) = [1,x]^T\) and \(q=1+p = 2\).
For the covariance function we choose \(\sigma^2c(\cdot,\cdot),\) where the correlation function \(c(\cdot,\cdot)\) has the Gaussian form described in the alternatives page on emulator prior correlation function (AltCorrelationFunction)
noting that in this case \(p=1\).
Estimation of the correlation length¶
We start with the estimation of the correlation length \(\delta\). In this example we will use the value of \(\delta\) that maximises the posterior distribution \(\pi^*_{\delta}(\delta)\), assuming that there is no prior information on \(\delta\), i.e. \(\pi(\delta)\propto \mathrm{const}\). The expression that needs to be maximised, according to the procedure page on building a Gaussian process emulator for the core problem (ProcBuildCoreGP), is
where
\(H\) is defined as
\(A\) is the \(n\times n\) correlation matrix of the design points, with elements \(c(x_i,x_j)\), \(i,j \in \{1,\cdots,n\}\).
Recall that in the above expressions the only term that is a function of \(\delta\) is the correlation matrix \(A\).
Reparameterisation¶
In order to maximise the posterior \(\pi^*_{\delta}(\delta)\), we reparameterise it using \(\tau = \ln(2\delta)\), to obtain \(\pi^*(\tau) = \pi^*_{\delta}(\exp(\tau/2))\). This has the benefit of making the optimisation problem unconstrained, because \(\delta\in (0,\infty)\), while \(\tau\in (-\infty,\infty)\).
Posterior function¶
Figure 2 shows \(\ln(\pi^*(\tau))\) as a function of \(\tau\).
Figure 2: The log posterior of \(\tau\) as a function of \(\tau\)
The maximum of this value can be obtained with any maximisation algorithm. In our case, we used Nelder - Mead. The value of \(\tau\) that maximises the posterior is -2.76 and the respective value of \(\delta\) is 0.25. From now on, we will refer to this value as \(\hat{\delta}\).
Because we have scaled the input to lie in \([0,1]\), \(\hat\delta\) is one quarter of the range of \(x\) over which we are fitting the emulator. This is a fairly typical value for a smoothness parameter, indicating that there is structure to the way that the output responds to this input in addition to the linear mean function that has been fitted, but that the response is not particularly “wiggly”. In terms of the original input scale, \(\hat\delta\) corresponds to a smoothness parameter of \(200\times 0.25 = 50\).
Estimates for the remaining parameters¶
Apart from the correlation lengths, the two other parameters of the Gaussian Process are \(\beta\) and \(\sigma^2\). Having estimated the correlation lengths, the estimate for \(\widehat\sigma^2\) is given by the equation above and
Note that in these equations, the matrix \(A\) is calculated using \(\hat{\delta}\). The values we got for the above parameters are \(\hat{\beta} = [ -47.30, 53.79]^T\) and \(\hat{\sigma}^2 = 92.89\) . Therefore, the fitted underlying linear trend is \(y=-47.3+53.79x\) (remembering that \(x\) is transformed to lie in [0,1]), with deviation from this line measured by a standard deviation of \(\sqrt {92.89}=9.64\).
Posterior mean and Covariance functions¶
The expressions for the posterior mean and covariance functions according to ProcBuildCoreGP are
and
Figure 3 shows the predictions of the emulator for 100 points uniformly spaced in \(1075, 1325\) in the original scale. The continuous line is the output of the simulator and the dashed line is the emulator’s output \(m^*\). The shaded areas represent 2 times the standard deviation of the emulator’s prediction, which is the square root of the diagonal of matrix \(v^*\).
Figure 3: Simulator (continuous line), emulator’s mean (dashed line) and 95% confidence intervals (shaded area)
We see that the true simulator output lies outside the two standard deviation interval for inputs around 1160 but otherwise seems to capture the behaviour of the simulator well. The next step would be to validate the emulator, at which point the anomaly around 1160 might be found, leading to a rebuilding and improvement of the emulator.
Validation¶
In this section we validate the above emulator according to the procedure page on validating a Gaussian process emulator (ProcValidateCoreGP).
The first step is to select the validation design. Following the advice in the discussion page on the design of a validation sample (DiscCoreValidationDesign), we chose \(n'=3p = 3\). We choose 2 points close to the original training points and one more halfway between 1140 and 1180, where the emulator makes the biggest jump. In the original input space of the simulator, the validation points we choose are
and in the transformed space
The simulator’s output at these points is
We then calculate the mean \(m^*(\cdot)\) and variance \(v^*(\cdot,\cdot)\) of the emulator at the validation design point \(D'\). The means are
and the standard deviations are
The predictions with error bars at 2 standard deviations are shown in Figure 4, which shows that all the predictions are within 2 standard deviations, with the exception of the prediction at \(\tilde{x}=1160\).
Figure 4: Emulator’s predictions at the validation points, with error bars at \(\pm 2\) standard deviations
Figure 5: Individual standardised errors for the prediction at the validation points
Figure 5, shows the individual standardised errors, indicating that they are all within the boundaries of \(\pm 2\) standard deviations, with the exception of the prediction at 1160, which is at 3 standard deviations.
The Mahalanobis distance diagnostic for this validation was found to be
when its theoretical mean is
The conclusion is that the emulator failed in both diagnostics, because the Mahalanobis distance is too large, and one of the individual standardised errors is close to 3. We will therefore rebuild the emulator using both the previous training points and the validation points for training.
Rebuilding the emulator¶
We now rebuild the emulator using both the training and the validation points. That is,
which in the transformed space are
The simulator’s output at these points is
Following the same procedure as before, we find the following values for the parameters:
Figure 6 shows the predictions of the emulator for 100 points uniformly spaced in 1075, 1325 in the original scale. A comparison with Figure 3 shows that the new emulator has a smaller predictive variance, and the mean of the prediction is closer to the actual value of the simulator.
Figure 6: Mean and 95% confidence intervals for the rebuilt emulator
Re-validation¶
We finally carry out the two validation tests for the new emulator. Using the same reasoning as before, we select the following validation points
which in the transformed input space are
The simulator’s output at these points is
The following the same new validation points we select are
The simulator’s output for these values is
The mean and standard deviation of the emulator at the new validation points is
Figure 7 shows the individual standardised errors. Note that they are now all within 2 standard deviations.
The Mahalanobis distance is now \(M=6.17\), whereas its theoretical mean is
and its variance is
Therefore, the Mahalanobis distance found is a typical value from its reference distribution, and taking also into account the fact that all the standardised errors are within 2 standard deviations, the emulator is declared valid.
Final Build¶
After the validation being successful, we can rebuild the emulator using all the data, which would provide us with its final version. The final set of 12 training points is
which in the transformed space are
The simulator’s output at these points is
The hyperparameter estimates obtained using the above training points are
Finally, figure 8 shows the predictions of the emulator for 100 points, uniformly spaced in 1075, 1325 in the original scale. Comparing with figures 3 and 6, we see that the addition of the extra training points has further reduced the predictive variance and moved the mean closer to the simulator’s output.
Figure 8: Mean and 95% confidence intervals for the final build of the emulator