sigmahatsquared | 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 |
A |
Correlation matrix (eg |
Ainv |
Inverse of the correlation matrix (eg |
d |
Vector of observations |
The formula is
\frac{y^T\left(A^{-1}-A^{-1}H(H^T A^{-1}H)^{-1} H^T
A^{-1}\right)y}{n-q-2}
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.
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)
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.