rtmvn {tmvmixnorm}

Random number generation for truncated multivariate normal distribution subject to linear inequality constraints

Description

rtmvn simulates truncated multivariate (p-dimensional) normal distribution subject to linear inequality constraints. The constraints should be written as a matrix (D) with lower and upper as the lower and upper bounds for those constraints respectively. Note that D can be non-full rank, which generalize many traditional methods.

Usage

rtmvn(
  n,
  Mean,
  Sigma,
  D = diag(1, length(Mean)),
  lower,
  upper,
  int = NULL,
  burn = 10,
  thin = 1
)

Arguments

n

number of random samples desired (sample size).

Mean

mean vector of the underlying multivariate normal distribution.

Sigma

positive definite covariance matrix of the underlying multivariate normal distribution.

D

matrix or vector of coefficients of linear inequality constraints.

lower

vector of lower bounds for truncation.

upper

vector of upper bounds for truncation.

int

initial value vector for Gibbs sampler (satisfying truncation), if NULL then determine automatically.

burn

burn-in iterations discarded (default as 10).

thin

thinning lag (default as 1).

Value

rtmvn returns a (n*p) matrix (or vector when n=1) containing random numbers which approximately follows truncated multivariate normal distribution.

Examples

# Example for full rank with strong dependence
d <- 3
rho <- 0.9
Sigma <- matrix(0, nrow=d, ncol=d)
Sigma <- rho^abs(row(Sigma) - col(Sigma))

D1 <- diag(1,d) # Full rank

set.seed(1203)
ans.1 <- rtmvn(n=1000, Mean=1:d, Sigma, D=D1, lower=rep(-1,d), upper=rep(1,d),
int=rep(0,d), burn=50)

apply(ans.1, 2, summary)

# Example for non-full rank
d <- 3
rho <- 0.5
Sigma <- matrix(0, nrow=d, ncol=d)
Sigma <- rho^abs(row(Sigma) - col(Sigma))

D2 <- matrix(c(1,1,1,0,1,0,1,0,1),ncol=d)
qr(D2)$rank # 2

set.seed(1228)
ans.2 <- rtmvn(n=100, Mean=1:d, Sigma, D=D2, lower=rep(-1,d), upper=rep(1,d), burn=10)

apply(ans.2, 2, summary)