# Sampling Random Numbers From The Truncated Multivariate Normal Distribution With Linear Constraints This function generates random numbers from the truncated multivariate normal distribution with mean equal to `mean` and covariance matrix `sigma` and general linear constraints [REMOVE_ME]$$ lower \le D x \le upperlower <= D x <= upper [REMOVE_ME_2]$$ with either rejection sampling or Gibbs sampling. ## 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 upperlower <= D x <= upper $$ with either rejection sampling or Gibbs sampling. ```r 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 - `n`: Number of random points to be sampled. Must be an integer $>= 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))`. - `D`: 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` ## Details This method allows for $r > d$ linear constraints, whereas `rtmvnorm` requires a full-rank matrix D $(d x d)$ and can only handle $r <= d$ constraints at the moment. The lower and upper bounds `lower` and `upper` are $(r x 1)$, the matrix `D` is $(r x d)$ and x is $(d x 1)$. The default case is $r = d$ and $D = I_d$. ## Warning This method will be merged with `rtmvnorm` in one of the next releases. ## Author(s) Stefan Wilhelm ## See Also `rtmvnorm` ## Examples ```r ## 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") ## End(Not run) ```

rtmvnorm2 function