# mvndst: Multivariate Normal Distribution CDF and Its Derivative In pedmod: Pedigree Models

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

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

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)){
function(par){
set.seed(1)
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_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.