lrstat (version 0.3.3)

pmvnormr: Multivariate Normal Probabilities over a Hyperrectangle

Description

Computes the probability that a multivariate normal random vector falls within a rectangular region defined by lower and upper bounds.

Usage

pmvnormr(
  lower = NULL,
  upper = NULL,
  mean = NULL,
  sigma,
  n0 = 1024,
  n_max = 16384,
  R = 8,
  abseps = 1e-04,
  releps = 0,
  seed = 314159,
  parallel = TRUE,
  nthreads = 0
)

Value

The estimated probability with the following attributes:

• method: "exact" for analytic methods or "qmc" for the Monte Carlo approach.

• error: The estimated error (half-width of 95% confidence interval).

• nsamples: The total number of samples used across all replications.

Arguments

lower

A numeric vector of lower integration limits.

upper

A numeric vector of upper integration limits.

mean

The mean vector. If NULL (default), a zero vector of appropriate length is used.

sigma

The covariance (or correlation) matrix of the distribution.

n0

Initial number of samples per replication for the Monte Carlo integration.

n_max

Maximum number of samples allowed per replication.

R

Number of independent replications used to estimate the error.

abseps

Absolute error tolerance for the probability calculation.

releps

Relative error tolerance for the probability calculation.

seed

Random seed for reproducibility. If 0, a seed is generated from the computer clock.

parallel

Logical; if TRUE, computations are performed in parallel.

nthreads

Number of threads for parallel execution. If 0, the default RcppParallel behavior is used.

Author

Kaifeng Lu,

Details

The function automatically selects the most efficient computation method based on the input:

• Analytic Methods: Used for univariate cases or multivariate distributions with a compound symmetry correlation structure and non-negative correlations.

• Monte Carlo Estimation: Used for all other cases. The algorithm employs a randomized quasi-Monte Carlo (QMC) approach using generalized Halton sequences.

To improve efficiency and accuracy, the QMC approach incorporates:

• Sequential Conditioning: Mimics the standardization and transformation approach used in mvtnorm::lpmvnorm, reducing the \(J\)-dimensional integral to a \((J-1)\)-dimensional problem over a hypercube.

• Adaptive Sampling: The number of samples per replication increases dynamically until the estimated error falls below abseps or releps, or until n_max is reached.

The standard error is derived from \(R\) independent replications. For high-dimensional problems, computations can be accelerated by setting parallel = TRUE, which distributes the replications across multiple CPU threads via nthreads.

Positive semidefinite sigma matrices (including singular cases) are supported in the general covariance branch via minimal diagonal stabilization during factorization.

Examples

# Example 1: Compound symmetry covariance structure and analytic method
n <- 5
mean <- rep(0, n)
lower <- rep(-1, n)
upper <- rep(3, n)
sigma <- matrix(0.5, n, n)
diag(sigma) <- 1
pmvnormr(lower, upper, mean, sigma)

# Example 2: General covariance structure and Monte Carlo method
n <- 5
mean <- rep(0, n)
lower <- rep(-1, n)
upper <- rep(3, n)
sigma <- matrix(c(1, 0.5, 0.3, 0.2, 0.1,
                  0.5, 1, 0.4, 0.3, 0.2,
                  0.3, 0.4, 1, 0.5, 0.3,
                  0.2, 0.3, 0.5, 1, 0.4,
                  0.1, 0.2, 0.3, 0.4, 1), nrow = n)
pmvnormr(lower, upper, mean, sigma, seed = 314159, nthreads = 1)