Multivariate normal integration

Overview

Many quantities reduce to the chance that a normal random vector falls inside a region: the best accuracy of a two-class classifier, or an error rate for a decision boundary 2012.14331 §1 2012.14331 §3.1. For a quadratic boundary there is a dedicated route. For any other shape the paper reports no closed form and no standard method or software 2012.14331 §1.

The idea all variants share: whiten the vector with the symmetric square root of its covariance, so the distribution becomes a standard normal and the region is mapped with it 2012.14331 §2.1. The integral is then taken over that standard normal.

Variants

Generalized chi-square (quadratic domains)

The region is written as a quadratic inequality in the standardised variable. An eigen-decomposition of its quadratic part turns the quadratic into a sum of scaled non-central chi-squared terms plus a normal term. The probability is then a one-dimensional integral of the generalized chi-square cdf at zero 2012.14331 §2.1. A flat boundary needs only a single normal cdf 2012.14331 §2.1. The paper names Ruben's, Imhof's, IFFT and ray methods for that last integral 2012.14331 §2.1. For one dimension and an interval, the paper writes the Bayes error with non-central chi-squared cdfs 2012.14331 §3.1.1. Card: Methods to integrate multinormals and compute classification measures.

Ray trace (arbitrary domains)

Standardise the domain, then look at every direction from the centre. Along one ray the domain is described by its sign at minus infinity and its crossing points 2012.14331 §2.2.2. Each direction contributes a weighted tail term built from the chi-squared survival function, and the probability is the integral of these terms over directions, divided by the total solid angle 2012.14331 §2.2.3 2012.14331 §2.2.4. Up to four dimensions the paper integrates over a grid on half the sphere. Above four dimensions it draws random uniform directions instead 2012.14331 §2.2.4. Complement, intersection and union of domains are applied to the ray traces directly 2012.14331 §2.2.5. For a function of a normal vector, the same traces give the cdf, and the pdf follows by numerical differentiation 2012.14331 §2.2.6.

Two limits to note. For a non-quadratic domain the sums run over crossings within a finite radius, and the mass beyond that radius is lost. The paper gives no numeric radius 2012.14331 §2.2.2. For classes defined as a minimum of pairwise quadratic regions, the paper says only the ray method computes the error rates 2012.14331 §3.4.

Monte Carlo (baseline)

Sampling points from the normal and counting how many fall in the region is the reference method in the paper's tests, not a method it proposes 2012.14331 §5. In a three-dimensional benchmark with equal covariances and 10^8 samples, the relative error in d' grows past about 10. Both integrators stay at machine precision up to about 75 2012.14331 §5 2012.14331 fig. 8a.

When to use which

The paper's results are author-reported, on the authors' MATLAB implementation and one hardware setup 2012.14331 §5. Re-measure before relying on the timings. For discriminability measures built on these integrals, see discriminability-index.

Implementation sketch

Common core of the variants. The quadratic branch and the ray branch differ in how they handle the domain. Dimension and sign conventions follow the paper.

# Probability that x ~ N(mu, Sigma) lies in a domain  [arxiv:2012.14331 §2]
# Input:  mu (k,), sigma (k, k), domain (quadratic, or a membership test)
# Output: p (scalar)

def normal_domain_probability(mu, sigma, domain):
    # Whitening matrix: symmetric square root of sigma  [arxiv:2012.14331 §2.1]
    sd = symmetric_sqrt(sigma)
    if domain.is_quadratic:
        # Standardise the quadratic, then eigen-decompose its quadratic part
        # [arxiv:2012.14331 §2.1]
        q2, q1, q0 = standardise_quadratic(domain, mu, sd)
        w, r = eigen(q2)
        # One 1-D integral of the generalized chi-square cdf at zero
        # [arxiv:2012.14331 §2.1]
        return gen_chi_square_prob_positive(w, r, q1, q0)
    # Otherwise average tail terms over directions  [arxiv:2012.14331 §2.2.4]
    directions, weights = half_sphere_quadrature(k)
    total = 0
    for n, weight in zip(directions, weights):
        # Sign at minus infinity and crossings along this ray
        # [arxiv:2012.14331 §2.2.2]
        sign0, roots = trace_ray(domain, mu, sd, n)
        total += weight * ray_alpha(sign0, roots, k)
    # Divide by the total solid angle  [arxiv:2012.14331 §2.2.3]
    return total / solid_angle(k)

Where the variants differ: the quadratic branch needs no sampling and no direction grid. The ray branch needs a root finder per direction and a finite search radius 2012.14331 §2.2.2. Monte Carlo replaces the whole function with the fraction of samples that land in the domain.

Contradictions

None known. Only one card covers this slug, so no two papers can be compared. The paper's text and its released code differ in several defaults, such as the dimension where the code switches from grid to Monte Carlo. These differences are listed on the card and are not resolved here.

Papers

• Methods to integrate multinormals and compute classification measures: the generalized chi-square and ray-trace methods for normal probabilities, and their comparison with Monte Carlo sampling.