stage1: Stage 1,2 and 3 optimization on toy dataset
In RobinHankin/calibrator: Bayesian Calibration of Complex Computer Codes

Description Usage Arguments Details Author(s) References See Also Examples

Perform O'Hagan's three stage optimization on the toy dataset. Function stage1() and stage2() find the optimal values for the hyperparameters and stage3() finds the optimal values for the three parameters.

stage1(D1, y, H1,  maxit,  trace=100, method="Nelder-Mead",
      directory = ".", do.filewrite=FALSE, do.print=TRUE,
      phi.fun, lognormally.distributed=FALSE, include.prior=TRUE, phi)
stage2(D1, D2, H1, H2, y, z, maxit, trace=100, method = "Nelder-Mead",
      directory = ".", do.filewrite=FALSE, do.print=TRUE,  extractor,
      phi.fun, E.theta, Edash.theta, isotropic=FALSE,
      lognormally.distributed = FALSE, include.prior = TRUE,
      use.standin = FALSE, rho.eq.1 = TRUE, phi) 
stage3(D1, D2, H1, H2, d, maxit, trace=100, method="Nelder-Mead",
      directory = ".", do.filewrite=FALSE, do.print=TRUE,
      include.prior = TRUE, lognormally.distributed=FALSE,
      theta.start=NULL, phi)

`maxit`	Maximum number of iterations as passed to `optim()`
`trace`	Amount of information displayed, as passed to `optim()`
`D1`	Matrix whose rows are points at which code output is known
`D2`	Matrix whose rows are points at which observations were made
`H1,H2`	Regressor basis functions for `D1` and `D2`
`y`	Code outputs. Toy example is `y.toy`
`z`	Observations. Toy example is `z.toy`
`d`	Data vector consisting of the code runs and observations
`extractor`	extractor function for `D1`
`E.theta,Edash.theta`	Expectation WRT theta, and dashed theta. Toy examples are `E.theta.toy()` and `Edash.theta.toy()`
`phi.fun`	Function to create hyperparameters; passed (in `stage1()` and `stage2()`) to `phi.change()`. Toy version is `phi.fun.toy()`
`method`	Method argument passed to `optim()`; qv
`include.prior`	Boolean variable with default `TRUE` meaning to include the prior distribution in the optimization process and `FALSE` meaning to use an uniformative prior (effectively uniform support). This variable is passed to `p.eqn4.supp()` for `stage1()`, `p.page4()` for `stage2()`, and `p.eqn8.supp()` for `stage3()`
`lognormally.distributed`	Boolean with `TRUE` meaning to use a lognormal distn. See `prob.theta` for details
`do.filewrite`	Boolean, with `TRUE` meaning to save a `load`able file `stage[123].<date>`, containing the interim value of `phi` and the corresponding optimand to `directory` at each evalution of the optimizer. If `FALSE`, don't write the files
`directory`	The directory to write files to; only matters if `do.filewrite` is `TRUE`
`isotropic`	In function `stage2()`, Boolean with default `FALSE` meaning to carry out a full optimization, and `TRUE` meaning to restrict the scope to isotroic roughness matrices. See details section below
`do.print`	Boolean, with default `TRUE` meaning to print interim values of `phi` at each evaluation
`use.standin`	In `stage2()`, a Boolean argument, with default `FALSE` meaning to use the real value for matrix `V.temp`, and `TRUE` meaning to use a standing that is the same size but contains fictitious values. The only time to set `use.standin` to `TRUE` is when debugging as it runs several orders of magnitude faster
`rho.eq.1`	Boolean, with default `TRUE` meaning to hold the value of `rho` constant at one (1)
`theta.start`	In `stage3()`, the starting point of the optimization with default `NULL` meaning to use the maximum likelihood point of the apriori distribution (ie `phi$theta.apriori$mean`)
`phi`	Hyperparameters. Used as initial values for the hyperparameters in the optimization routines

The three functions documented here carry out the multi-stage optimization detailed in KOH2001 (actually, KOH2001 only defined stage 1 and stage 2, which estimated the hyperparameters. What is here called “stage3()” estimates the true value of theta given the hyperparameters).

stage1() carries out stage 1 of KOH2001 which is used to estimate psi1 using optimization.

In function stage2(), setting argument isotropic to TRUE will force phi$omegastar_x to be a function of a length one scalar. The value of phi$omegastar_x used will depend on pdm.maker.psi2() (an internal function appearing in hpa.fun.toy()). In stage2(), several kludges are made. The initial conditions are provided by argument phi. The relevant part of this is phi$psi2.

Function stage2() estimates psi2 and rho and lambda, using optimization. Note that psi2 includes sigma2squared in addition to omegastar_X (in the toy case, psi2 has three elements: the first two are the diagonal of omegastar_x and the third is sigma2squared but this information is encoded in phi.fun.toy(), which changes from application to application).

Function stage3() attempts to find the maximum likelihood estimate of theta, given hyperparameters and observations, using optimization

Robin K. S. Hankin

M. C. Kennedy and A. O'Hagan 2001. Bayesian calibration of computer models. Journal of the Royal Statistical Society B, 63(3) pp425-464
M. C. Kennedy and A. O'Hagan 2001. Supplementary details on Bayesian calibration of computer models, Internal report, University of Sheffield. Available at http://www.tonyohagan.co.uk/academic/ps/calsup.ps
R. K. S. Hankin 2005. Introducing BACCO, an R bundle for Bayesian analysis of computer code output, Journal of Statistical Software, 14(16)

toys, phi.fun.toy

data(toys)
stage1(D1=D1.toy,y=y.toy,H1=H1.toy, maxit=5, phi.fun=phi.fun.toy, phi=phi.toy)

##now try with a slightly bigger dataset:
##Examples below take a few minutes to run:

set.seed(0)
data(toys)
jj <- create.new.toy.datasets(D1.toy , D2.toy)
y.toy <- jj$y.toy
z.toy <- jj$z.toy
d.toy <- jj$d.toy

phi.toy.stage1 <- stage1(D1=D1.toy, y=y.toy, H1=H1.toy, maxit=10, phi.fun=phi.fun.toy, phi=phi.toy)

phi.toy.stage2 <- stage2(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy,
 y=y.toy, z=z.toy, extractor=extractor.toy,
phi.fun=phi.fun.toy, E.theta=E.theta.toy, Edash.theta=Edash.theta.toy,
maxit=3, phi=phi.toy.stage1)

stage3(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, d=d.toy, maxit=3, phi=phi.toy.stage2)

# Now try with the true values of the hyperparameters:
phi.true <- phi.true.toy(phi=phi.toy)

stage3(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, d=d.toy, maxit=3, phi=phi.true)