mcCovMat: Monte Carlo estimation of a covariance matrix.

Description

Calculate an estimate of the covariance matrix for the parameter estimates of a db or beta binomial distribution via simulation.

Usage

mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)
# S3 method for mleDb
mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)
# S3 method for mleBb
mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)
# S3 method for Dbdpars
mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)
# S3 method for Bbdpars
mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)
# S3 method for default
mcCovMat(object, nsim = 100, seed=NULL, maxit=1000)

Arguments

object

An object of class either "mleDb", "mleBb", Dbdpars or Bbdpars. In the first two cases such an object would be returned by the function mleDb() or by mleBb(). In the second two cases such an object would be returned by the function makeDbdpars() or by makeBbdpars().

nsim

Integer scalar. The number of simulations to be used to produce the Monte Carlo estimate of the covariance matrix.

seed

Integer scalar. The seed for the random number generator. If not specified it is randomly sampled from the sequence 1:1e5.

maxit

Integer scalar. The maximum number of iterations to be undertaken by optim() when fitting models to the simulated data.

Value

A two-by-two positive definite (with any luck!) numeric matrix. It is an estimate of the covariance matrix of the parameter estimates.

It has an attribute "seed" which is the seed that was used for the random number generator. This is either the value of the argument seed or (if this argument was left NULL) the value that was randomly sampled from 1:1e5.

Details

The procedure is to simulate nsim data sets, all of the same size. This will be the size of the data set to which object was fitted), in the case of the "mleDb" and "mleBb" methods, and will be the value of the ndata argument supplied to the “make” function in the case of the "Dbdpars" and "Bbdpars" methods. The simulations are from models determined by the parameter value contained in object.

From each such simulated data, parameter estimates are obtained. The covariance matrix of these latter parameter estimates (adjusted for the fact that the true parameters are known in a simulation) is taken to be the required covariance matrix estimated.

The default method simply throws an error.

Examples

Run this code

# NOT RUN {
X   <- hmm.discnp::SydColDisc
X$y <- as.numeric(X$y)
X   <- split(X,f=with(X,interaction(locn,depth)))
x   <- X[[19]]$y
fit <- mleDb(x, ntop=5)
set.seed(42)
CM.m <- mcCovMat(fit,nsim=500) # Lots of simulations!
CM.a <- vcov(fit)
CM.n <- solve(nHess(fit,x))
cat("Monte Carlo:\n\n")
print(CM.m)
cat("Analytic:\n\n")
print(CM.a)
cat("Numeric:\n\n")
print(CM.n)
X     <- hrsRcePred
top1e <- X[X$sbjType=="Expert","top1"]
fit   <- mleBb(top1e,size=10)
CM.m  <- mcCovMat(fit,nsim=500) # Lots of simulations!
CM.a  <- vcov(fit)
CM.n  <- solve(nHess(fit,top1e))
cat("Monte Carlo:\n\n")
print(CM.m)
cat("Analytic:\n\n")
print(CM.a)
cat("Numeric:\n\n")
print(CM.n)
# }