Estimator for sigma squared

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

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)

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.