rtmvnorm2 | R Documentation |
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 <= D x <= upper
with either rejection sampling or Gibbs sampling.
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"), ...)
n |
Number of random points to be sampled. Must be an integer >= 1. |
mean |
Mean vector (d x 1), default is |
sigma |
Covariance matrix (d x d), default is |
lower |
Vector of lower truncation points (r x 1),
default is |
upper |
Vector of upper truncation points (r x 1),
default is |
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 |
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.
This method will be merged with rtmvnorm
in one of the next releases.
Stefan Wilhelm
rtmvnorm
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.