# sigmahatsquared: Estimator for sigma squared In emulator: Bayesian emulation of computer programs

## Description

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.

## Usage

 ```1 2``` ```sigmahatsquared(H, Ainv, d) sigmahatsquared.A(H, A, d) ```

## Arguments

 `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

## Details

The formula is

ommitted; see pdf

where y is the data vector, H the matrix whose rows are the regressor functions of the design matrix, A the correlation matrix, n the number of observations and q the number of elements in the basis function.

## Author(s)

Robin K. S. Hankin

## References

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

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30``` ```## 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) 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) ```

emulator documentation built on May 30, 2017, 12:06 a.m.