tmvnorm | R Documentation |
Density, distribution function and random generation for the multivariate truncated normal distribution
with mean vector mu
, covariance matrix sigma
, lower truncation limit lb
and upper truncation limit ub
.
The truncation limits can include infinite values. The Monte Carlo (type = "mc"
) uses a sample of size B
, while the
quasi Monte Carlo (type = "qmc"
) uses a pointset of size ceiling(n/12)
and estimates the relative error using 12 independent randomized QMC estimators.
n |
number of observations |
x , q |
vector of quantiles |
B |
number of replications for the (quasi)-Monte Carlo scheme |
log |
logical; if |
mu |
vector of location parameters |
sigma |
covariance matrix |
lb |
vector of lower truncation limits |
ub |
vector of upper truncation limits |
check |
logical; if |
... |
additional arguments, currently ignored |
type |
string, either of |
dtmvnorm
gives the density, ptmvnorm
and pmvnorm
give the distribution function of respectively the truncated and multivariate Gaussian distribution and rtmvnorm
generate random deviates.
dtmvnorm(x, mu, sigma, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4) ptmvnorm(q, mu, sigma, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4) rtmvnorm(n, mu, sigma, lb, ub)
Zdravko I. Botev, Leo Belzile (wrappers)
Z. I. Botev (2017), The normal law under linear restrictions: simulation and estimation via minimax tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.
d <- 4; lb <- rep(0, d)
mu <- runif(d)
sigma <- matrix(0.5, d, d) + diag(0.5, d)
samp <- rtmvnorm(n = 10, mu = mu, sigma = sigma, lb = lb)
loglik <- dtmvnorm(samp, mu = mu, sigma = sigma, lb = lb, log = TRUE)
cdf <- ptmvnorm(samp, mu = mu, sigma = sigma, lb = lb, log = TRUE, type = "q")
# Exact Bayesian Posterior Simulation Example
# Vignette, example 5
## Not run:
data("lupus"); # load lupus data
Y <- lupus[,1]; # response data
X <- as.matrix(lupus[,-1]) # construct design matrix
n <- nrow(X)
d <- ncol(X)
X <- diag(2*Y-1) %*% X; # incorporate response into design matrix
nusq <- 10000; # prior scale parameter
C <- solve(diag(d)/nusq + crossprod(X))
sigma <- diag(n) + nusq*tcrossprod(X) # this is covariance of Z given beta
est <- pmvnorm(sigma = sigma, lb = 0)
# estimate acceptance probability of crude Monte Carlo
print(attributes(est)$upbnd/est[1])
# reciprocal of acceptance probability
Z <- rtmvnorm(sigma = sigma, n = 1e3, lb = rep(0, n))
# sample exactly from auxiliary distribution
beta <- rtmvnorm(n = nrow(Z), sigma = C) + Z %*% X %*% C
# simulate beta given Z and plot boxplots of marginals
boxplot(beta, ylab = expression(beta))
# output the posterior means
colMeans(beta)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.