BACCO (version 1.0-50)

betahat.fun: Calculates MLE coefficients of linear fit

Description

Determines the MLE (least squares) regression coeffients for the specified regression basis.

The .A form needs only A (and not Ainv), thus removing the need to calculate a matrix inverse. Note that this form is slower than the other if Ainv is known in advance, as solve(.,.) is slow.

If Ainv is not known in advance, the two forms seem to perform similarly in the cases considered here and in the goldstein package.

Usage

betahat.fun(xold, Ainv, d, give.variance=FALSE, func)
betahat.fun.A(xold, A, d, give.variance=FALSE, func)

Arguments

xold

Data frame, each line being the parameters of one run

A

Correlation matrix, typically provided by corr.matrix()

Ainv

Inverse of the correlation matrix A

d

Vector of results at the points specified in xold

give.variance

Boolean, with TRUE meaning to return information on the variance of $\hat{\beta}$ and default FALSE meaning to return just the least squares estimator

func

Function to generate regression basis; defaults to regressor.basis

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

data(toy)
val <- toy
H <- regressor.multi(val)
d <- apply(H,1,function(x){sum((0:6)*x)})


fish <- rep(2,6)
A <- corr.matrix(val,scales=fish, power=2)
Ainv <- solve(A)

# now add suitably correlated Gaussian noise:
d <-  as.vector(rmvnorm(n=1,mean=d, 0.1*A))

betahat.fun(val , Ainv , d)           # should be close to c(0,1:6)


# Now look at the variances:
betahat.fun(val,Ainv,give.variance=TRUE, d)


     # now find the value of the a priori expectation (ie the regression
     # plane) at an unknown point:
x.unknown <- rep(0.5 , 6)
regressor.basis(x.unknown) %*% betahat.fun(val, Ainv, d)

     # compare the a-priori with the a-postiori:
interpolant(x.unknown, d, val, Ainv,scales=fish)
     # Heh, it's the same!  (of course it is, there is no error here!)


     # OK, put some error on the old observations:
d.noisy <- as.vector(rmvnorm(n=1,mean=d,0.1*A))

     # now compute the regression point:
regressor.basis(x.unknown) %*% betahat.fun(val, Ainv, d.noisy)

     # and compare with the output of interpolant():
interpolant(x.unknown, d.noisy, val, Ainv, scales=fish)
     # there is a difference!



     # now try a basis function that has superfluous degrees of freedom.
     # we need a bigger dataset.  Try 100:
val <- latin.hypercube(100,6)
colnames(val) <- letters[1:6]
d <- apply(val,1,function(x){sum((1:6)*x)})
A <- corr.matrix(val,scales=rep(1,6), power=1.5)
Ainv <- solve(A)

    
betahat.fun(val, Ainv, d, func=function(x){c(1,x,x^2)})
     # should be c(0:6 ,rep(0,6).  The zeroes should be zero exactly
     # because the original function didn't include any squares.


## And finally a sanity check:
betahat.fun(val, Ainv, d, func=function(x){c(1,x,x^2)})-betahat.fun.A(val, A, d, func=function(x){c(1,x,x^2)})