mvndst | R Documentation |
Provides an approximation of the multivariate normal distribution CDF over a hyperrectangle and the derivative with respect to the mean vector and the covariance matrix.
mvndst( lower, upper, mu, sigma, maxvls = 25000L, abs_eps = 0.001, rel_eps = 0L, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, method = 0L, n_sequences = 8L, use_tilting = FALSE ) mvndst_grad( lower, upper, mu, sigma, maxvls = 25000L, abs_eps = 0.001, rel_eps = 0L, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, method = 0L, n_sequences = 8L, use_tilting = FALSE )
lower |
numeric vector with lower bounds. |
upper |
numeric vector with upper bounds. |
mu |
numeric vector with means. |
sigma |
covariance matrix. |
maxvls |
maximum number of samples in the approximation. |
abs_eps |
absolute convergence threshold. |
rel_eps |
relative convergence threshold. |
minvls |
minimum number of samples. Negative values provides a default which depends on the dimension of the integration. |
do_reorder |
|
use_aprx |
|
method |
integer with the method to use. Zero yields randomized Korobov lattice rules while one yields scrambled Sobol sequences. |
n_sequences |
number of randomized quasi-Monte Carlo sequences to use. More samples yields a better estimate of the error but a worse approximation. Eight is used in the original Fortran code. If one is used then the error will be set to zero because it cannot be estimated. |
use_tilting |
|
mvndst:
An approximation of the CDF. The "n_it"
attribute shows the number of
integrand evaluations, the "inform"
attribute is zero if the
requested precision is achieved, and the "abserr"
attribute
shows 3.5 times the estimated standard error.
mvndst_grad:
A list with
likelihood
: the likelihood approximation.
d_mu
: the derivative with respect to the the mean vector.
d_sigma
: the derivative with respect to the covariance matrix
ignoring the symmetry (i.e. working the n^2 parameters with
n being the dimension rather than the n(n + 1) / 2
free parameters).
# simulate covariance matrix and the upper bound set.seed(1) n <- 10L S <- drop(rWishart(1L, 2 * n, diag(n) / 2 / n)) u <- rnorm(n) system.time(pedmod_res <- mvndst( lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n), maxvls = 1e6, abs_eps = 0, rel_eps = 1e-4, use_aprx = TRUE)) pedmod_res # compare with mvtnorm if(require(mvtnorm)){ mvtnorm_time <- system.time(mvtnorm_res <- mvtnorm::pmvnorm( upper = u, sigma = S, algorithm = GenzBretz( maxpts = 1e6, abseps = 0, releps = 1e-4))) cat("mvtnorm_res:\n") print(mvtnorm_res) cat("mvtnorm_time:\n") print(mvtnorm_time) } # with titling system.time(pedmod_res <- mvndst( lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n), maxvls = 1e6, abs_eps = 0, rel_eps = 1e-4, use_tilting = TRUE)) pedmod_res # compare with TruncatedNormal if(require(TruncatedNormal)){ TruncatedNormal_time <- system.time( TruncatedNormal_res <- TruncatedNormal::pmvnorm( lb = rep(-Inf, n), ub = u, sigma = S, B = attr(pedmod_res, "n_it"), type = "qmc")) cat("TruncatedNormal_res:\n") print(TruncatedNormal_res) cat("TruncatedNormal_time:\n") print(TruncatedNormal_time) } # check the gradient system.time(pedmod_res <- mvndst_grad( lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n), maxvls = 1e5, minvls = 1e5, abs_eps = 0, rel_eps = 1e-4, use_aprx = TRUE)) pedmod_res # compare with numerical differentiation. Should give the same up to Monte # Carlo and finite difference error if(require(numDeriv)){ num_res <- grad( function(par){ set.seed(1) mu <- head(par, n) S[upper.tri(S, TRUE)] <- tail(par, -n) S[lower.tri(S)] <- t(S)[lower.tri(S)] mvndst( lower = rep(-Inf, n), upper = u, sigma = S, mu = mu, maxvls = 1e4, minvls = 1e4, abs_eps = 0, rel_eps = 1e-4, use_aprx = TRUE) }, c(numeric(n), S[upper.tri(S, TRUE)]), method.args = list(d = .01, r = 2)) d_mu <- head(num_res, n) d_sigma <- matrix(0, n, n) d_sigma[upper.tri(d_sigma, TRUE)] <- tail(num_res, -n) d_sigma[upper.tri(d_sigma)] <- d_sigma[upper.tri(d_sigma)] / 2 d_sigma[lower.tri(d_sigma)] <- t(d_sigma)[lower.tri(d_sigma)] cat("numerical derivatives\n") print(rbind(numDeriv = d_mu, pedmod = pedmod_res$d_mu)) print(d_sigma) cat("\nd_sigma from pedmod\n") print(pedmod_res$d_sigma) # for comparison }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.