mvndst: Multivariate Normal Distribution CDF and Its Derivative

In pedmod: Pedigree Models

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	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 doi: 10.1111/rssb.12162.

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).

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
}

pedmod documentation built on Sept. 11, 2022, 5:05 p.m.