tmvtnorm (version 0.8-3)

rtmvnorm: Sampling Random Numbers From Truncated Multivariate Normal Distribution

Description

This function generates random numbers from the truncated multivariate normal 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

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

Arguments

n

Number of observations.

mean

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

sigma

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

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), the Fortan Gibbs sampler ("gibbs") and the old Gibbs sampler implementation in R ("gibbsR").

...

additional parameters for Gibbs sampling, given to rmvtnorm.gibbs, such as burn.in.samples and start.value

Details

The generation of random numbers from a truncated multivariate normal distribution is done using either rejection sampling or Gibbs sampler. Rejection sampling is done from the standard multivariate normal distribution. So we use the function rmvnorm of the mvtnorm package. In order to speed up the generation of N samples from the truncated distribution, we first calculate the acceptance rate alpha from the truncation points and then generate N/alpha samples iteratively until we have got N samples. This typically does not take more then 2-3 iterations. Rejection sampling may be very inefficient when the support region is small (i.e. in higher dimensions) which results in very low acceptance rates alpha. In this case the Gibbs sampler is preferable. The Gibbs sampler samples from univariate conditional distributions, so all samples can be accepted except for a burn-in period. The number of burn-in samples to be discarded can be specified, as well as a start value. If no start value is given we determine a start value from the support region using either lower bound or upper bound if they are finite, or 0 otherwise. The Gibbs sampler has been reimplemented in Fortran 90 for performance reasons (algorithm="gibbs"). The old R implementation is still accessible through algorithm="gibbsR".

References

Johnson, N./Kotz, S. (1970). Distributions in Statistics: Continuous Multivariate Distributions Wiley & Sons, pp. 70--73 Horrace, W. (2005). Some Results on the Multivariate Truncated Normal Distribution. Journal of Multivariate Analysis, 94, 209--221 Jayesh H. Kotecha and Petar M. Djuric (1999). Gibbs Sampling Approach For Generation of Truncated Multivariate Gaussian Random Variables IEEE Computer Society, 1757--1760

See Also

ptmvnorm, pmvnorm, rmvnorm, dmvnorm

Examples

dtmvnorm(x=c(0,0))
dtmvnorm(x=c(0,0), mean=c(1,1), upper=c(0,0))

###########################################
#
# Example 1: 
# rejection sampling        
#
############################################

sigma <- matrix(c(4,2,2,3), ncol=2)
x <- rtmvnorm(n=500, mean=c(1,2), sigma=sigma, upper=c(1,0))
plot(x, main="samples from truncated bivariate normal distribution",
  xlim=c(-6,6), ylim=c(-6,6), 
  xlab=expression(x[1]), ylab=expression(x[2]))
abline(v=1, lty=3, lwd=2, col="gray")
abline(h=0, lty=3, lwd=2, col="gray")

###########################################
#
# Example 2: 
# Gibbs sampler for 4 dimensions
#
############################################

C = matrix(0.8, 4, 4)
diag(C)=rep(1, 4)
lower = rep(-4, 4)
upper = rep(-1, 4)

# acceptance rate alpha
alpha = pmvnorm(lower=lower, upper=upper, mean=rep(0,4), sigma=C)
alpha

# Gibbs sampler
X1=rtmvnorm(n=20000, mean = rep(0,4), sigma=C, lower=lower, upper=upper, 
  algorithm="gibbs", burn.in.samples=100)
# Rejection sampling
X2=rtmvnorm(n=5000, mean = rep(0,4), sigma=C, lower=lower, upper=upper)

colMeans(X1)
colMeans(X2)

plot(density(X1[,1]), col="red", lwd=2, main="Gibbs vs. Rejection")
lines(density(X2[,1]), col="blue", lwd=2)
legend("topleft",legend=c("Gibbs Sampling","Rejection Sampling"), 
  col=c("red","blue"), lwd=2)