rtmvnorm2: Sampling Random Numbers From The Truncated Multivariate Normal Distribution With Linear Constraints

Description

This function generates random numbers from the truncated multivariate normal distribution with mean equal to mean and covariance matrix sigma and general linear constraints $$lower \le D x \le upper$$ with either rejection sampling or Gibbs sampling.

Usage

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

Arguments

Number of random points to be sampled. Must be an integer $\ge 1$.

mean

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

sigma

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

lower

Vector of lower truncation points (r x 1), default is rep( Inf, length = length(mean)).

upper

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

Matrix for linear constraints (r x d), defaults to diagonal matrix (d x d), i.e. r = d.

algorithm

Method used, possible methods are the Fortan Gibbs sampler ("gibbs", default), the Gibbs sampler implementation in R ("gibbsR") and rejection sampling ("rejection")

…

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

Warning

This method will be merged with rtmvnorm in one of the next releases.

Details

This method allows for $r > d$ linear constraints, whereas rtmvnorm requires a full-rank matrix D $(d \times d)$ and can only handle $r \le d$ constraints at the moment. The lower and upper bounds lower and upper are $(r \times 1)$, the matrix D is $(r \times d)$ and x is $(d \times 1)$. The default case is $r = d$ and $D = I_d$.

Examples

Run this code

# NOT RUN {
################################################################################
#
# Example 5a: Number of linear constraints r > dimension d
#
################################################################################

# general linear restrictions a <= Dx <= b with x (d x 1); D (r x d); a,b (r x 1)

# Dimension d=2, r=3 linear constraints
#
# a1 <=    x1 + x2 <= b2
# a2 <=    x1 - x2 <= b2
# a3 <= 0.5x1 - x2 <= b3
#
# [ a1 ] <= [ 1     1 ] [ x1 ] <= [b1]
# [ a2 ]    [ 1    -1 ] [ x2 ]    [b2]
# [ a3 ]    [ 0.5  -1 ]           [b3]

D <- matrix(
      c(  1,  1,
          1, -1,
        0.5, -1), 3, 2, byrow=TRUE)
a <- c(0, 0, 0)
b <- c(1, 1, 1)

# mark linear constraints as lines
plot(NA, xlim=c(-0.5, 1.5), ylim=c(-1,1))
for (i in 1:3) {
  abline(a=a[i]/D[i, 2], b=-D[i,1]/D[i, 2], col="red")
  abline(a=b[i]/D[i, 2], b=-D[i,1]/D[i, 2], col="red")
}

### Gibbs sampling for general linear constraints a <= Dx <= b
mean <- c(0, 0)
sigma <- matrix(c(1.0, 0.2, 
                  0.2, 1.0), 2, 2)
x0 <- c(0.5, 0.2) # Gibbs sampler start value                  
X <- rtmvnorm2(n=1000, mean, sigma, lower=a, upper=b, D, start.value=x0)

# show random points within simplex
points(X, pch=20, col="black")
# }