mvt.mcmc: Estimate Parameters of a Multivariate t Distribution Using the MCMC

Description

Use the MCMC to obtain estimate of parameters of a multivariate t distribution.

Usage

mvt.mcmc(X, prior.lower.v, prior.upper.v,
           prior.Mu0=rep(0, ncol(X)), prior.Sigma0=diag(10000, ncol(X)),
           prior.p=ncol(X), prior.V=diag(1, ncol(X)),
           initial.v=NULL, initial.Sigma=NULL,
           nmcmc=10000, nburn=nmcmc/10, nthin=1, seed=1)

Arguments

a matrix of observations with one subject per row.

prior.lower.v

lower bound of degrees of freedom (df).

prior.upper.v

upper bound of df.

prior.Mu0

mean vector of multivariate normal prior of the location. The default value is 0.

prior.Sigma0

variance matrix of multivariate normal prior of the location. The default value is a diagonal matrix with diagonal entries equal to 10000.

prior.p

the df of wishart prior of inverse of the scale matrix. The default value is dimensions of observation.

prior.V

the scale matrix of wishart prior of inverse of the scale matrix. The default value is identity matrix.

initial.v

the initial value of the df. The default is NULL. For the default, the value will be generated by using the ECME Algorithm.

initial.Sigma

the initial value of the scale matrix. The default is NULL. For the default, the value will be generated by using the ECME Algorithm.

nmcmc

number of iterations. The default value is 10000.

nburn

number of burn-in. The default value is nmcmc/10.

nthin

output every nthin-th sample. The default value is 1 (no thinning).

seed

random seed. The default value is 1. Note that seed can be an integer from 1 to 14. The internal state of the OpenBUGS random number generator can be set to one of 14 predefined states.

Value

Mu.save

a matrix of locations of the distribution, one row per iteration.

Sigma.save

a three dimensional array of scale matrix of the distribution. Sigma.save[,,i] is the result from the ith iteration.

v.save

a vector of df of the distribution, one component per iteration.

Details

To generate samples from the full conditional distribution of df, the slice sampling was used and the code was adapted from http://www.cs.toronto.edu/~radford/ftp/slice-R-prog.

Examples

Run this code

# NOT RUN {
  mu1 <- mu2 <- sigma12 <- sigma22 <- 100
  rho12 <- 0.9
  Sigma <- matrix(c(sigma12, rho12*sqrt(sigma12*sigma22),
                    rho12*sqrt(sigma12*sigma22), sigma22),
                  nrow=2)
  k <- 8
  N <- 100
  X <- rmvt(N, sigma=Sigma, df=k, delta=c(mu1, mu2))

  result <- mvt.mcmc(X, 4, 25)
  colMeans(result$Mu.save)
  apply(result$Sigma.save, c(1,2), mean)
  mean(result$v.save)
# }