sigmahatsquared {emulator} | R Documentation |
Returns maximum likelihood estimate for sigma squared. The
“.A
” form does not need Ainv
, thus removing the need to
invert A
. Note that this form is slower than
the other if Ainv
is known in advance, as solve(.,.)
is slow.
sigmahatsquared(H, Ainv, d) sigmahatsquared.A(H, A, d)
H |
Regressor matrix (eg as returned by regressor.multi() ) |
A |
Correlation matrix (eg corr.matrix(val) ) |
Ainv |
Inverse of the correlation matrix (eg solve(corr.matrix(val)) ) |
d |
Vector of observations |
Robin K. S. Hankin
J. Oakley and A. O'Hagan, 2002. “Bayesian Inference for the Uncertainty Distribution of Computer Model Outputs”, Biometrika 89(4), pp769-784
R. K. S. Hankin 2005. “Introducing BACCO, an R bundle for Bayesian analysis of computer code output”, Journal of Statistical Software, 14(16)
## First, set sigmasquared to a value that we will try to estimate at the end: REAL.SIGMASQ <- 0.3 ## First, some data: val <- latin.hypercube(100,6) H <- regressor.multi(val,func=regressor.basis) ## now some scales: fish <- c(1,1,1,1,1,4) ## A and Ainv A <- corr.matrix(as.matrix(val),scales=fish,power=2) Ainv <- solve(A) ## a real relation; as used in helppage for interpolant: real.relation <- function(x){sum( (1:6)*x )} ## use the real relation: d <- apply(val,1,real.relation) ## now add some Gaussian process noise: d.noisy <- as.vector(rmvnorm(n=1,mean=d, REAL.SIGMASQ*A)) ## now estimate REAL.SIGMASQ: sigmahatsquared(H,Ainv,d.noisy) ## That shouldn't be too far from the real value specified above. ## Finally, a sanity check: sigmahatsquared(H,Ainv,d.noisy) - sigmahatsquared.A(H,A=A,d.noisy)