V.fun {calibrator}

Variance matrix for observations

Description

Determines the variance/covariance matrix for the observations and code run points.

Usage

V.fun(D1, D2, H1, H2,  extractor,
E.theta, Edash.theta, give.answers=FALSE, test.for.symmetry=FALSE, phi)

Arguments

D1

Matrix of code run points

D2

Matrix of observation points

H1

Regression function for D1

H2

Regression function for D2

extractor

Function to extract x.star and t.vec from D1

Edash.theta

Function to return expectation of H with respect to the alternative distribution of \theta; Edash.theta.toy is the example for the toy dataset

E.theta

Expectation of h WRT theta over the apriori distribution. Note that this function must be updated if h1() changes.

give.answers

Boolean (defaulting to FALSE) with TRUE meaning to return a list whose elements are V and its constituent parts, viz line1 to line6. This argument is used mainly for debugging.

test.for.symmetry

Boolean with TRUE meaning to calculate each element of C explicitly, and default FALSE meaning to calculate only the elements of C that lie on or over the diagonal and use the fact that C is symmetric to calculate the other matrix elements. For n observations, this means n(n+1)/2 evaluations, compared with n^2 for the full case. The time saving is considerable, even for small matrices.

Set this argument to TRUE only when debugging, or testing accuracy

phi

Hyperparameters

Details

See KOH2001 for full details on page 3 of the supplement

Value

If give.answers is the default value of FALSE, returns a matrix of covariances for use in p.page4().

If give.answers is TRUE, returns a named list of (currently) 17 elements. Elements one to six are lines one to six respectively from page 3 of the supplement; subsequent lines give intermediate steps in the calculation. The final element is the matrix is the covariances as returned when give.answers is FALSE.

Note

This function takes a long time to run

Author(s)

Robin K. S. Hankin

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

See Also

tt.fun,p.page4

Examples

data(toys)
(jj <-V.fun(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, 
  extractor=extractor.toy, 
  Edash.theta=Edash.theta.toy,
  E.theta=E.theta.toy,  phi=phi.toy))


## Now note that V.fun() changes with the PRIOR used for theta:
phi.different.theta <-  phi.change(old.phi=phi.toy,
     theta.apriori.mean=c(100,100,100),phi.fun=phi.fun.toy)
V.fun(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, 
  extractor=extractor.toy, 
  Edash.theta=Edash.theta.toy,
  E.theta=E.theta.toy,  phi=phi.different.theta)
## different!


## Now compare jj above with V.fun() calculated with
## different phi2:

phi.toy.new <- phi.change(phi.fun=phi.fun.toy, old.phi = phi.toy, psi2=c(8,8,8))

V.fun(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, 
  extractor=extractor.toy, 
  Edash.theta=Edash.theta.toy,
  E.theta=E.theta.toy,  phi=phi.toy.new)

## different!



## Not run: 
data(toys)
set.seed(0)
jj <- create.new.toy.datasets(D1=D1.toy , D2=D2.toy)
y.toy <- jj$y.toy
z.toy <- jj$z.toy
d.toy <- jj$d.toy

v.fun <- function(...){V.fun(D1=D1.toy, D2=D2.toy, H1=H1.toy, H2=H2.toy, 
     extractor=extractor.toy, Edash.theta=Edash.theta.toy,
     E.theta=E.theta.toy, phi=phi.toy, give=TRUE)}

Rprof(file="~/f.txt");ignore <- v.fun();Rprof(file=NULL)
system("cd ; R CMD Rprof ~/f.txt > ~/ff.txt")

## End(Not run)