Estimate Parameters of a Multivariate t Distribution Using the MCMC

Share:

Description

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

Usage

1
2
3
4
  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

X

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.

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.

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.

See Also

mvt.ecme

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
## 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, 10000, 4, 25)
  colMeans(result$Mu.save)
  apply(result$Sigma.save, c(1,2), mean)
  mean(result$v.save)

## End(Not run)