rtmvt {tmvtnorm}

Sampling Random Numbers From The Truncated Multivariate Student t Distribution

Description

This function generates random numbers from the truncated multivariate Student-t distribution with mean equal to mean and covariance matrix sigma, lower and upper truncation points lower and upper with either rejection sampling or Gibbs sampling.

Usage

rtmvt(n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)), df = 1, 
  lower = rep(-Inf, length = length(mean)), 
  upper = rep(Inf, length = length(mean)),
  algorithm=c("rejection", "gibbs"), ...)

Arguments

n

Number of random points to be sampled. Must be an integer >= 1.

mean

Mean vector, default is rep(0, length = ncol(x)).

sigma

Covariance matrix, default is diag(ncol(x)).

df

Degrees of freedom parameter (positive, may be non-integer)

lower

Vector of lower truncation points,\ default is rep(-Inf, length = length(mean)).

upper

Vector of upper truncation points,\ default is rep( Inf, length = length(mean)).

algorithm

Method used, possible methods are rejection sampling ("rejection", default) and the R Gibbs sampler ("gibbs").

...

additional parameters for Gibbs sampling, given to the internal method rtmvt.gibbs(), such as burn.in.samples, start.value and thinning, see details

Details

We sample x \sim T(\mu, \Sigma, df) subject to the rectangular truncation lower \le x \le upper. Currently, two random number generation methods are implemented: rejection sampling and the Gibbs Sampler.

For rejection sampling algorithm="rejection", we sample from rmvt and retain only samples inside the support region. The acceptance probability will be calculated with pmvt. pmvt does only accept integer degrees of freedom df. For non-integer df, algorithm="rejection" will throw an error, so please use algorithm="gibbs" instead.

The arguments to be passed along with algorithm="gibbs" are:

burn.in.samples

number of samples in Gibbs sampling to be discarded as burn-in phase, must be non-negative.

start.value

Start value (vector of length length(mean)) for the MCMC chain. If one is specified, it must lie inside the support region (lower \le start.value \le upper). If none is specified, the start value is taken componentwise as the finite lower or upper boundaries respectively, or zero if both boundaries are infinite. Defaults to NULL.

thinning

Thinning factor for reducing autocorrelation of random points in Gibbs sampling. Must be an integer \ge 1. We create a Markov chain of length (n*thinning) and take only those samples j=1:(n*thinning) where j %% thinning == 0 Defaults to 1 (no thinning of the chain).

Warning

The same warnings for the Gibbs sampler apply as for the method rtmvnorm.

Author(s)

Stefan Wilhelm <[e-mail]>, Manjunath B G <[e-mail]>

References

Geweke, John F. (1991) Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints. Computer Science and Statistics. Proceedings of the 23rd Symposium on the Interface. Seattle Washington, April 21-24, 1991, pp. 571–578 An earlier version of this paper is available at https://www.researchgate.net/publication/2335219_Efficient_Simulation_from_the_Multivariate_Normal_and_Student-t_Distributions_Subject_to_Linear_Constraints_and_the_Evaluation_of_Constraint_Probabilities

Examples

###########################################################
#
# Example 1
#
###########################################################	

# Draw from multi-t distribution without truncation
X1 <- rtmvt(n=10000, mean=rep(0, 2), df=2)
X2 <- rtmvt(n=10000, mean=rep(0, 2), df=2, lower=c(-1,-1), upper=c(1,1))

###########################################################
#
# Example 2
#
###########################################################	

df = 2
mu = c(1,1,1)
sigma = matrix(c(  1, 0.5, 0.5,
                 0.5,   1, 0.5,
                 0.5, 0.5,   1), 3, 3)
lower = c(-2,-2,-2)
upper = c(2, 2, 2)

# Rejection sampling
X1 <- rtmvt(n=10000, mu, sigma, df, lower, upper)

# Gibbs sampling without thinning
X2 <- rtmvt(n=10000, mu, sigma, df, lower, upper, 
  algorithm="gibbs")

# Gibbs sampling with thinning
X3 <- rtmvt(n=10000, mu, sigma, df, lower, upper, 
  algorithm="gibbs", thinning=2)	
   
plot(density(X1[,1], from=lower[1], to=upper[1]), col="red", lwd=2,
  main="Gibbs vs. Rejection")
lines(density(X2[,1], from=lower[1], to=upper[1]), col="blue", lwd=2)
legend("topleft",legend=c("Rejection Sampling","Gibbs Sampling"), 
  col=c("red","blue"), lwd=2)

acf(X1)  # no autocorrelation in Rejection sampling
acf(X2)  # strong autocorrelation of Gibbs samples
acf(X3)  # reduced autocorrelation of Gibbs samples after thinning