# stage1

##### Stage 1,2 and 3 optimization on toy dataset

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.

- Keywords
- array

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

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

##### 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-464M. 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.psR. K. S. Hankin 2005.

*Introducing BACCO, an R bundle for Bayesian analysis of computer code output*, Journal of Statistical Software, 14(16)

##### See Also

##### Examples

```
# NOT RUN {
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)
# }
```

*Documentation reproduced from package calibrator, version 1.2-8, License: GPL-2*