
sample.from.exp.est()
samples from the
appropriate distribution. This function is not designed for the general user: it is tailored for
use in the environment of the National Oceanographic Centre, with a
particular version of the specialist model
makeinputfiles(number.of.runs = 100, gaussian = TRUE,
directoryname="~/goldstein/genie-cgoldstein/", filename="QWERTYgoin", expert.estimates, area.outside=0.05)
sample.from.exp.est(number.of.runs, expert.estimates, gaussian=TRUE, area.outside=0.05)
upper
and lower
columns are used: here these values data(expert.estimates)
to load
a sample dataset that was supplied by Bob Marsha
and b
as the $2.5%\rm{ile}$
and $97.5%\rm{ile}$ respectively./working/jrd/sat/rksh/goldstein
. The database
results.table
is made using the shell scripts currently in
/users/sat/rksh/goldstein/emulator
. Note that makeinputfiles(number.of.runs=n)
creates files
numbered from $0$ to $n-1$: so be careful of
off-by-one errors. It's probably best to avoid reference to the
The suffix number of a file matches the name of its tmp
file
(so, for example, file with suffix number 15 writes output to files
tmp/tmp.15
and tmp/tmp.avg.15
).
expert.estimates
, results.table
data(expert.estimates) system("mkdir /users/sat/rksh/tmp",ignore=TRUE)
makeinputfiles(number.of.runs = 100, gaussian = TRUE,
directoryname="~/tmp/", expert.estimate=expert.estimates)
data(results.table)
data(expert.estimates)
output.col <- 25
wanted.row <- 1:27
wanted.cols <- c(2:9,12:19)
val <- results.table[wanted.row , wanted.cols]
mins <- expert.estimates$low
maxes <- expert.estimates$high
normalize <- function(x){(x-mins)/(maxes-mins)}
unnormalize <- function(x){mins + (maxes-mins)*x}
jj <- t(apply(val,1,normalize))
jj <- as.data.frame(jj)
names(jj) <- names(val)
val <- jj
scales.optim <- exp(c( -2.63, -3.03, -2.24, 2.61,
-1.65, -3.13, -3.52, 3.16, -3.32, -2.53, -0.25, -2.55, -4.98, -1.59,
-4.40, -0.81))
d <- results.table[wanted.row , output.col]
A <- corr.matrix(val, scales=scales.optim, power=1.5)
Ainv <- solve(A)
x <- sample.from.exp.est(1000,exp=expert.estimates)
x <- t(apply(x,1,normalize))
ensemble <- interpolant.quick(x , d , val , Ainv, scales=scales.optim, power=1.5)
hist(ensemble)
Run the code above in your browser using DataLab