test.csv: To cross-validate a Gaussian Process emulator
In stilt: Separable Gaussian Process Interpolation (Emulation)

Description Usage Arguments Details Value References See Also Examples

To test an emulator using cross-validation. Any reasonable number of runs can be excluded from the emulator for testing. Then, the emulator is used to predict at the excluded parameter settings, and the results can be plotted. Optionally, prediction standard deviation can also be plotted. The program is reproducible if the seed number is specified.

1 2	test.csv(final.emul, num.test, plot.std, theseed = NULL, test.runind = NULL, make.plot = TRUE)

`final.emul`	A standard emulator object. It can be generated by, for example, the 'emulator' function.
`num.test`	Number of test runs to withhold at random from the emulator. Can not be more than half of the ensemble.
`plot.std`	If `TRUE`, prediction std. is plotted around the emulator prediction
`theseed`	Seed for random number generator. Default is `NULL`.
`test.runind`	If not `NULL`, a monotonically increasing vector of run indices at which to test the emulator. Has to contain `num.test` number of elements. If `NULL`, run indices are generated at random. Default is `NULL`.
`make.plot`	If `TRUE` (the default) produced a cross-validation plot.

The function withholds model runs from the emulator, and then uses the withheld emulator to predict at the excluded parameter settings. Any reasonable number of runs up to 1/2 of the ensemble can be excluded. If make.plot=TRUE the plot of model output and emulator prediction at the excluded setting is produced. Setting plot.std=TRUE adds the prediction standard deviation to the plot. If plot.std=TRUE but make.plot=FALSE a warning is given and nothing is plotted. The default setting is to exclude num.test runs at random. If theseed is specified the results become reproducible as the seed passed on to the R random number generator is fixed. Optionally, the excluded run indices can be specified explicitly via the vector test.runind. If both theseed and test.runind arguments are specified, then theseed is ignored.

The function prints out the excluded run indices and outputs some basic guidance. The function handles runs which are at parameter space boundaries by passing NA to the relevant rows of output matrices. If prediction at a particular withheld run gives error (most likely due to the fact that extrapolation is not allowed and the excluded run is at the corner of parameter space) that particular run is not plotted. Model output is beige, and emulator output is brown.

By 'withholding' runs from an emulator it is meant that the emulator statistical parameters are kept the same, whereas the withheld runs are eliminated from all relevant emulator components. The ensemble size final.emul$p is reduced accordingly.

List with components

`$model.out.test`	Model output at test parameter settings [row,col] = [test run index, time index]
`$emul.out.test`	Emulator output at test parameter settings [row,col] = [test run index, time index]
`$emul.std.test`	Emulator prediction standard deviation at test parameter settings [row,col] = [test run index, time index]
`$coverage`	Empirical 95% coverage, expressed as a fraction

The time vector for columns of all matrix components is final.emul$t.vec

R. Olson and W. Chang (2013): Mathematical framework for a separable Gaussian Process Emulator. Tech. Rep., available from
www.scrimhub.org/resources/stilt/Olson_and_Chang_2013_Stilt_Emulator_Technical_Report.pdf.

test.all, predict.emul

# Test the 1D emulator by withholding the second parameter
# setting, and then predicting at the withheld run
data(emul.1D)
test.1D <- test.csv(final.emul=emul.1D, num.test=1, plot.std=FALSE,
                    test.runind=2, make.plot=FALSE)

# Create a custom plot
plot.default(NA, xlim=range(emul.1D$t.vec), ylim=range(test.1D$model.out.test),
             xlab="Year", ylab="Sample Output", main="Emulator Test at Theta=1",
             cex.lab=1.3, cex.axis=1.3, cex.main=1.3)
lines(emul.1D$t.vec, test.1D$model.out.test, lwd=4, lty=2, col="green")
# The emulator prediction is close to perfect! Yay!
lines(emul.1D$t.vec, test.1D$emul.out.test, lwd=2, col="yellow")

# Produce a legend
model.max <- max(test.1D$model.out.test)
# We already know that time ranges from 0 to 10 so figuring out x placement
# for the legend is not hard
legend(1, model.max, c("Model Output", "Emulator Predictions"),
       col=c("green", "yellow"), lty=c(2,1), lwd=c(4,2), cex=1.3)



# Test SICOPOLIS emulator at three pre-selected parameter settings
# specified by theseed=1234321
data(emul.Sicopolis)
# Plot the standard deviation in all three cases
test.csv(final.emul=emul.Sicopolis, num.test=3, plot.std=TRUE, theseed=1234321,
         make.plot=TRUE)