Last chance! 50% off unlimited learning
Sale ends in
var.conditional() calculates the conditional variance-covariance
matrix, and cond.sample() samples from the appropriate
multivariate t distribution.cond.sample(n = 1, x, xold, d, A, Ainv, scales = NULL, pos.def.matrix =
NULL, func = regressor.basis, ...)
var.conditional(x, xold, d, A, Ainv, scales = NULL, pos.def.matrix = NULL,
func = regressor.basis, distance.function = corr, ...)cond.sample(), the number of observations
to take, defaulting to 1corr())corr()cond.sample() returns a $n\times p$ matrix
whose rows are independent samples from the appropriate multivariate
$t$ distribution. Here, $p$ is the number of rows of x
(ie the number of simulated design points). Consider a case where there
are just two simulated design points, close to each other but far from
any point of the original design points. Then function
cond.sample(n=4, ...) will give four numbers which are close to
one another but have high (between-instantiation) variance.Function var.conditional() calculates the denominator of equation
3 of Oakley and OHagan 2002. This function is intended to be called by
cond.sample() but might be interesting per se when
debugging or comparing different choices of simulated design points.
x; uncertainty in x is captured by assuming it to
be drawn from a pdf X.
The basic idea is to estimate $m^*$ at simulated design
points using cond.sample(), which samples from the multivariate
t distribution conditional on the data d at the design points.
The random datavector of estimates $m^*$ is called ddash. We repeat this process many times, each time estimating
$\eta(\cdot)$ using the augmented dataset
c(d,ddash) as a training set.
For each estimated $\eta(\cdot)$, we have a complete emulator that can be used to build up an ensemble of estimates.
R. K. S. Hankin 2005.
regressor.basis, for a more visually informative
example of cond.sample() et seq; and interpolant
for more examples# Now we use the functions. First we set up some design points:
# Suppose we are given the toy dataset and want to know the PDF of
# fourth power of the response at point x, where uncertainty in x
# may be represented as it being drawn from a normnl distribution
# with mean c(0.5,0.5,...,0.5) and a variance of 0.001.
data(toy)
val <- toy
real.relation <- function(x){sum( (0:6)*x )}
H <- regressor.multi(val)
d <- apply(H,1,real.relation)
# and some scales (which are assumed to be known):
fish <- rep(1,6)
fish[6] <- 4
# And determine A and Ainv:
A <- corr.matrix(val,scales=fish, power=2)
Ainv <- solve(A)
# and add some suitably correlated Gaussian noise:
d.noisy <- as.vector(rmvnorm(n=1, mean=d, 0.1*A))
names(d.noisy) <- names(d)
# Now some simulation design points. Choose n'=6:
xdash <- matrix(runif(36),ncol=6)
xdash <- rbind(xdash,xdash[6,])
colnames(xdash) <- colnames(val)
# And just for fun, we insert a useless seventh simulation design point
# (it is useless because it is a copy of the sixth). We do this
# in order to test the coding:
rownames(xdash) <- c("alpha","beta","gamma","delta","epsilon","zeta","zeta.copy")
# Print the variance matrix:
var.conditional(x=xdash,xold=val,d=d.noisy,A=A,Ainv=Ainv,scales=fish)
# Note that the sixth and seventh columns are identical (and so,
# therefore, are the sixth and seventh rows) as expected.
# Now sample from the conditional t distribution. Taking n=3 samples:
cond.sample(n=3, x=xdash, xold=val, d=d.noisy, A=A, Ainv=Ainv, scales = fish,
func = regressor.basis)
# Note the last two columns are identical, as expected.
# Just as a test, what is the variance matrix at the design points?
var.conditional(x=val,xold=val,d=d.noisy,A=A,Ainv=Ainv,scales=fish)
# (this should be exactly zero)
# Next, we apply the methods of OO2002 using Monte Carlo techniques.
# We will generate 10 different versions of eta:
number.of.eta <- 10
# And, for each eta, we will sample from the posterior t distribution 11 times:
number.of.X <- 11
# create an augmented design matrix, of the design points plus the
# simulated design points:
design.augmented <- rbind(val,xdash)
Ainv.augmented <- corr.matrix(design.augmented, scales=fish)
out <- NULL
for(i in 1:number.of.eta){
# Create random data by sampling from the conditional multivariate t
# at the simulated design points xdash, from the t-distribution given the data d:
ddash <- cond.sample(n=1, x=xdash, xold=val, d=d.noisy, Ainv=Ainv, scales=fish)
# Now use the emulator to calculate m^* at points chosen from the
# PDF of X:
jj <-
interpolant.quick(x=rmvnorm(n=number.of.X,rep(0.5,6),diag(6)/1000),
d=c(d.noisy,ddash),
xold=design.augmented,
Ainv=Ainv.augmented,
scales=fish)
out <- c(out,jj)
}
# histogram of the fourth power:
hist(out^4, col="gray")
# See oo2002 for another example of cond.sample() in useRun the code above in your browser using DataLab