Distribution Function for a Normal Order Statistic

Computes lower- or upper-tail probabilities for the \(r\)th order statistic of a multivariate normal random vector.

pordmvnormr(
  z,
  r,
  mean = NULL,
  sigma,
  lower.tail = TRUE,
  n0 = 1024,
  n_max = 16384,
  R = 8,
  abseps = 1e-04,
  releps = 0,
  seed = 314159,
  parallel = TRUE,
  nthreads = 0
)

Arguments

z: The threshold value.
r: The order statistic index, with \(1 \le r \le m\).
mean: The mean vector. If NULL (default), a zero vector of appropriate length is used.
sigma: The covariance (or correlation) matrix of the distribution.
lower.tail: Logical; if TRUE (default), probabilities are \(P(Z_{(r)} \le z)\), otherwise \(P(Z_{(r)} \ge z)\).
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.

Value

The estimated probability with the following attributes:

method: "analytic" for analytic or one-dimensional integration methods, or "qmc" if randomized quasi-Monte Carlo was used for at least one inclusion-exclusion term.
error: The accumulated estimated error from the multivariate normal probability calculations.
nsamples: The total number of samples used across all randomized quasi-Monte Carlo terms.

Details

Let \(Z_{(1)} \le \cdots \le Z_{(m)}\) denote the ordered components of \(Z\). The lower-tail probability is evaluated as \(P(Z_{(r)} \le z) = P(\sum_i 1\{Z_i \le z\} \ge r)\). The upper-tail probability is evaluated as \(P(Z_{(r)} \ge z) = P(\sum_i 1\{Z_i \ge z\} \ge m-r+1)\).

For compound-symmetry covariance matrices with non-negative correlation, the function uses a one-dimensional conditioning integral. If the component means are equal, the conditional tail is an ordinary binomial tail; otherwise it is a Poisson-binomial tail.

For general covariance matrices, the function uses the equivalent inclusion-exclusion representation over upper-orthant multivariate normal probabilities. This reduces each rectangle probability term to a lower dimensional subvector probability.

Author

Kaifeng Lu, kaifenglu@gmail.com

Examples


n <- 5
mean <- rep(0, n)
sigma <- matrix(0.5, n, n)
diag(sigma) <- 1
pordmvnormr(z = 1, r = 4, mean = mean, sigma = sigma, nthreads = 1)
#> [1] 0.7891928
#> attr(,"method")
#> [1] "analytic"
#> attr(,"error")
#> [1] 0
#> attr(,"nsamples")
#> [1] 1
pordmvnormr(z = 1, r = 4, mean = mean, sigma = sigma,
            lower.tail = FALSE, nthreads = 1)
#> [1] 0.2108072
#> attr(,"method")
#> [1] "analytic"
#> attr(,"error")
#> [1] 0
#> attr(,"nsamples")
#> [1] 1