In RobinHankin/emulator: Bayesian Emulation of Computer Programs

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

Overview

To cite the emulator package in publications please use Hankin 2005. The emulator package provides R-centric functionality for working with Gaussian processes. The focus is on approximate evaluation of complex computer codes. The package is part of the the BACCO suite of software.

Installation

You can install the released version of emulator from CRAN with:

# install.packages("emulator")  # uncomment this to use the package
library("emulator")

The package is maintained on github.

The `emulator` package in use

set.seed(0)
options(digits=3)
fish <- c(1,1,4)  # roughness scales
var <- 0.3      # variance
real.relation <- function(x){sum( (1:3)*x )} # real relation
val  <- latin.hypercube(7,3)
colnames(val) <- c("alpha","beta","gamma")
A <- corr.matrix(val,scales=fish)
d <-  as.vector(rmvnorm(n=1,mean=apply(val,1,real.relation),var*A))

Suppose we have a complicated computer program which takes three parameters as input, and we can run it a total of seven times at different points in parameter space:

val
d

Above, val shows the seven points in parameter space at which we have run the code, and d shows the output at those points. Now suppose we wish to know what the code would have produced at point $p=(0.5, 0.5, 0.5)$, at which the point has not actually been run. This is straightforward with the package:

p <- c(0.5,0.5,0.)
fish <- c(1,1,4)
A <- corr.matrix(val,scales=fish)
interpolant(p, d, val, A = A, scales=fish, give=TRUE)

Above, object fish is a vector of roughness length ("scales") corresponding to the small-scale covariance properties of our function. This may be estimated from the problem or from the datapoints. Matrix A is a normalized variance-covariance matrix for the points of val.

The output gives various aspects of the Gaussian process associated with the original observations. The most interesting one is mstar.star which indicates that the best estimate for the code's output, if it were to be run at point $p$, would be about 2.41.