BACCO (version 1.0-50)

stage1: Stage 1,2 and 3 optimization on toy dataset

Description

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.

Usage

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, 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)

Arguments

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

Regressor basis functions for D1

H2

Regressor basis functions for 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

Expectation WRT theta. Toy example is E.theta.toy()

Edash.theta

Expectation WRT the dashed theta. Toy example is 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

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 loadable file stage[123]., containing the interim value of phi and the corresponding optimand to directory at each evalution of the opt

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.

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

Details

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 $\psi_1$ 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 $\psi_2$ and $\rho$ and $\lambda$, using optimization. Note that $\psi_2$ includes $\sigma_2^2$ in addition to omegastar_X (in the toy case, $\psi_2$ has three elements: the first two are the diagonal of omegastar_x and the third is $\sigma_2^2$ 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

References

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.shef.ac.uk/~st1ao/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)

See Also

toys, phi.fun.toy

Examples

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)