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, niter, 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)

Arguments

a matrix of observations with one subject per row.

niter

number of iterations.

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.

Value

Mu.savea matrix of locations of the distribution, one row per iteration.
Sigma.savea three dimensional array of scale matrix of the distribution. Sigma.save[,,i] is the result from the ith iteration.
v.savea 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

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, 10000, 4, 25)
  colMeans(result$Mu.save)
  apply(result$Sigma.save, c(1,2), mean)
  mean(result$v.save)