rtmvn: Random number generation for truncated multivariate normal...

In tmvmixnorm: Sampling from Truncated Multivariate Normal and t Distributions

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)

tmvmixnorm documentation built on Sept. 19, 2020, 1:07 a.m.