Function to create (phi.fun.toy) or modify
(phi.change) toy hyperparameters \(\phi\)
in a form suitable for passing to the other functions in the library.
The user should never make \(\phi\) by hand; always use one of these functions
phi.fun.toy(rho, lambda, psi1, psi1.apriori, psi2, psi2.apriori,
theta.apriori)
phi.change(phi.fun, old.phi = NULL, rho = NULL, lambda = NULL,
psi1 = NULL, psi1.apriori=NULL, psi1.apriori.mean=NULL,
psi1.apriori.sigma=NULL, psi2 = NULL, psi2.apriori=NULL,
psi2.apriori.mean=NULL, psi2.apriori.sigma=NULL,
theta.apriori=NULL, theta.apriori.mean=NULL,
theta.apriori.sigma=NULL)In phi.change(), the name of the function that
creates the hyperparameters. Use phi.fun.toy() for the toy
dataset
In function phi.change(), the hyperparameter
object \(\phi\) to be modified
Correlation hyperparameter appearing in main equation
Noise hyperparameter
Roughness lengths hyperparameter for design matrix
D1. Internal function pdm.maker.psi1() takes
psi1 as an argument and returns omega_x,
omega_t and sigma1squared.
Recall that \(\Omega_x\) and \(Omega_t\) are
arbitrary functions of \(\psi_1\). In this case,
the values are omega_x=psi1[1:2], omega_t=psi1[3:4]
and sigma1squared=psi1[6]
A priori PDF for \(\psi_1\). In the form
of a two element list with first element (mean) the mean and
second element (sigma) the covariance matrix; distribution of
the logarithms is assumed to be multivariate normal. In the toy
example, the mean is a vector of length six (the first five are
\(\psi_1\) and the sixth is for
\(\sigma_1^2\)), and the variance is the
corresponding six-by-six matrix. Use function prob.psi1() to
calculate the apriori probability density for a particular value of
\(\psi_1\)
In function phi.change.toy(), use this
argument to change just the mean of psi1 (and leave the value
of sigma unchanged)
In function phi.change.toy(), use
this argument to change just the variance matrix of psi1
Roughness lengths hyperparameter for D2.
Internal function pdm.maker.psi2() takes psi2 as an
argument and returns omegastar_x and sigma2squared.
In phi.fun.toy(), the values are omegastar_x=psi2[1:2]
and sigma2squared=psi2[3].
NB: function stage2() optimizes more than just psi2.
It simultaneously optimizes psi2 and lambda and
rho
A priori PDF for \(\psi_2\) and hyperparameters \(\rho\) and \(\lambda\) (in that order).
As for psi1.apriori, this is in the form of a list with the
first element (mean) the mean and second element
(sigma) the covariance matrix; the logs are multivariate
normal. In the toy example, the mean is a vector of length five.
The first and second elements of the mean are the apriori mean of
\(\rho\) and \(\lambda\) respectively; the third
and fourth elements are the apriori mean of \(\psi_2\) (that
is, x and y respectively); and the fifth is the mean
of \(\sigma_2^2\).
The second element of phi.toy$psi2.apriori, sigma, is
the corresponding four-by-four variance matrix. Use function
prob.psi2() to calculate the apriori probability density of a
particular value of \(\psi_2\)
In phi.change.toy(), use to change
just the mean of psi2
In
phi.change.toy(), use to change just the variance matrix of
psi2
Apriori PDF for
\(\theta\). As above, in the form of a list with elements
for the mean and covariance. The distribution is multivariate
normal (NB: The distribution is multivariate normal and NOT
lognormal! To be explicit: \(\log(\theta)\) is
lognormally distributed). Use function prob.theta() to
calculate the apriori probability density of a particular value of
\(\theta\)
In phi.change.toy(), use to change
just the mean of theta
In phi.change.toy(),
use to change just the variance matrix of theta
Returns a list of several elements:
Correlation hyperparameter
Noise hyperparameter
Roughness lengths hyperparameter for D1
Apriori mean and variance matrix for psi1
Roughness lengths hyperparameter for D2
Apriori mean and variance matrix for psi2
Apriori mean and variance matrix for the parameters
Positive definite matrix for the lat/long part of
D1, whose diagonal is psi1[1:2]
Positive definite matrix for the code parameters theta,
whose diagonal is psi1[3:5]
Positive definite matrix for use in equation 13 of
the supplement; represents distances between rows of D2
variance
variance
Upper triangular Cholesky decomposition for omega_x
Lower triangular Cholesky decomposition for omega_x
Upper triangular Cholesky decomposition for omega_t
Lower triangular Cholesky decomposition for omega_t
Precalculated matrix for use in
Edash.theta(...,fast.but.opaque=TRUE)
Precalculated matrix for use in
Edash.theta(...,fast.but.opaque=TRUE)
Precalculated scalar for use in
ht.fun(...,fast.but.opaque=TRUE)
Precalculated scalarfor use in
tt.fun()
Upper triangular Cholesky decomposition for A
Lower triangular Cholesky decomposition for A
Note that this toy function contains within itself
pdm.maker.toy() which extracts omega_x and
omega_t and sigma1squared from psi1.
This will need to be changed for real-world applications.
Earlier versions of the package had pdm.maker.toy()
defined separately.
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)
# NOT RUN {
phi.fun.toy(100,101,1:6,list(mean=rep(1,6),sigma=1+diag(6)),50:55,
list(mean=rep(0,4),sigma=0.1+diag(4)),
list(mean=0.1+(1:3),sigma=2.1+diag(3)))
phi.fun.toy(rho=1, lambda=1,
psi1 = structure(c(1.1, 1.2, 1.3, 1.4, 1.5, 0.7),
.Names = c("x", "y", "A","B", "C","s1sq")),
psi1.apriori = list(
mean=rep(0,6), sigma=0.4+diag(6)),
psi2=structure(c(2.1, 2.2), .Names = c("x","y")),
psi2.apriori = list(mean=rep(0,5),sigma=0.2+diag(5)),
theta.apriori = list(mean=0.1+(1:3),sigma=2.1+diag(3))
)
data(toys)
phi.change(phi.fun=phi.fun.toy, old.phi = phi.toy, rho = 100)
phi.change(phi.fun=phi.fun.toy, old.phi = phi.toy,
theta.apriori.sigma = 4*diag(3))
identical(phi.toy, phi.change(phi.fun=phi.fun.toy, old.phi=phi.toy))
# }
Run the code above in your browser using DataCamp Workspace