---
kind: concept
slug: multivariate-normal-integration
title: Multivariate normal integration
summary: >-
  Computes the probability that a normal random vector lies in a region. Use a
  generalized chi-square route when the boundary is quadratic, and a ray-trace
  integral over directions for any other shape. Monte Carlo sampling is the
  simple baseline, but it loses accuracy once the probability is tiny.
tags: [probability, numerical-integration, classification]
papers: ["2012.14331"]
related: [discriminability-index]
status: draft
updated_at: 2026-10-19
---

## 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 [arxiv:2012.14331 §1] [arxiv: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 [arxiv: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 [arxiv: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
[arxiv:2012.14331 §2.1]. A flat boundary needs only a single normal cdf
[arxiv:2012.14331 §2.1]. The paper names Ruben's, Imhof's, IFFT and ray methods
for that last integral [arxiv:2012.14331 §2.1]. For one dimension and an
interval, the paper writes the Bayes error with non-central chi-squared cdfs
[arxiv:2012.14331 §3.1.1]. Card: [Methods to integrate multinormals and compute classification measures](/papers/2012.14331).

### 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 [arxiv: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
[arxiv:2012.14331 §2.2.3] [arxiv: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 [arxiv:2012.14331 §2.2.4]. Complement,
intersection and union of domains are applied to the ray traces directly
[arxiv: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
[arxiv: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 [arxiv: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 [arxiv: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
[arxiv: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 [arxiv:2012.14331 §5]
[arxiv:2012.14331 fig. 8a].

## When to use which

| Situation | Use | Why |
|---|---|---|
| Region bounded by a quadratic, such as the Bayes boundary of two normals | Generalized chi-square | Reduces to one 1-D integral after an eigen-decomposition [arxiv:2012.14331 §2.1]. The paper reports it as the fastest of its methods in its benchmark [arxiv:2012.14331 §5] |
| Any other shape, or classes defined by several quadratic regions | Ray trace | No closed form needed. For minimum-of-quadratics classes, the paper says only the ray method computes the error rates [arxiv:2012.14331 §3.4] |
| Dimension up to four, and reproducible results wanted | Ray trace on a grid | Deterministic quadrature over half the sphere [arxiv:2012.14331 §2.2.4] |
| Dimension above four | Ray trace with random directions | The paper switches to random directions, so the estimate is stochastic [arxiv:2012.14331 §2.2.4] |
| Quick estimate, probabilities that are not tiny, no implementation time | Monte Carlo | Simple to write. In some of the paper's cases it ran faster than the ray trace, so runtime depends on the problem [arxiv:2012.14331 §5] [arxiv:2012.14331 fig. 8b] |
| Very small probabilities, or d' above about 10 | Not Monte Carlo | Sampling loses accuracy there [arxiv:2012.14331 §5] [arxiv:2012.14331 fig. 8a] |

The paper's results are author-reported, on the authors' MATLAB implementation
and one hardware setup [arxiv:2012.14331 §5]. Re-measure before relying on the
timings. For discriminability measures built on these integrals, see
[discriminability-index](/concepts/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.

```pseudocode
# 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 [arxiv: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](/papers/2012.14331) and are not resolved here.

## Papers

- [Methods to integrate multinormals and compute classification measures](/papers/2012.14331): the generalized chi-square and ray-trace methods for normal probabilities, and their comparison with Monte Carlo sampling.