pedmod (version 0.2.4)

mvndst: 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
)

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 the \(n^2\) parameters with \(n\) being the dimension rather than the \(n(n + 1) / 2\) free parameters).

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

TRUE if a heuristic variable reordering should be used. TRUE is likely the best value.

use_aprx

TRUE if a less precise approximation of pnorm and qnorm should be used. This may reduce the computation time while not affecting the result much.

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

TRUE if the minimax tilting method suggested by Botev (2017) should be used. See tools:::Rd_expr_doi("10.1111/rssb.12162").

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
}