mvndst {pedmod} | R Documentation |
Multivariate Normal Distribution CDF and Its Derivative
Description
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.
Usage
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
)
Arguments
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 |
|
Value
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 then^2
parameters withn
being the dimension rather than then(n + 1) / 2
free parameters).
Examples
# 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
}