Methods to integrate multinormals and compute classification measures

In plain words

Many models of decisions assume that the measured quantities follow a bell-shaped (normal) distribution. To check how well such a model does, you need the chance that an outcome lands in some region, for example "the score is above this curved boundary". For most shapes there is no formula for that, and counting points on a grid or drawing random samples is slow or inaccurate when the chance is very small.

The authors give two numerical routes. For regions cut out by a quadratic curve or surface, they turn the problem into a sum of known distributions. For any other region, they send rays out from the centre of the distribution, note where each ray enters and leaves the region, and add up the probability along each ray over all directions. These tools give the best possible accuracy of a two-class classifier, and they let you compare that with simpler indices.

The software is a Matlab toolbox. Its results agree with random sampling where sampling works, and they stay accurate in cases where sampling breaks down.

TL;DR

For a normal vector, the probability of any domain can be computed from a generalized chi-square form when the boundary is quadratic, and by a ray-trace integral over angles when it is not §2.1 §2.2. The same integrals give Bayes-optimal error rates, a discriminability index d'_b, and checks of whether a normal model fits sample data §3. Author-reported: the toolbox agrees with Monte Carlo where sampling works and stays accurate to d' about 75, where sampling fails beyond about 10 §5 fig. 8a. The code is MIT-licensed LICENSE.

Problem

Normal models are used throughout vision, neuroscience, robotics, psychology and economics, and their predicted performance needs the normal integrated over specific domains §1. The best accuracy of a classifier with a quadratic boundary is such an integral §1 §3.1. No general closed form, standard numerical method or software existed for arbitrary normals and domains §1. Grids waste effort on regions of low density or outside the domain. Monte Carlo cannot reach a chosen precision directly and is unreliable when the probability is tiny, as for highly separated classes §1.

Equal-variance shortcuts are common, but they limit how models can be tested and compared §1. The paper also addresses the discriminability indices that are used for unequal variances, which differ from the Bayes optimum §3.1.4.

Method

The paper works in standardised coordinates. Let \bm{x} \sim N(\bm{\mu}, \bm{\Sigma}) and \mathbf{S} = \bm{\Sigma}^{1/2} be the symmetric square root, the "sd matrix" §2.1. A quadratic domain is

q(\bm{x}) = \bm{x}' \mathbf{Q}_2 \bm{x} + \bm{q}_1' \bm{x} + q_0 > 0 \qquad [\text{eq. 1}]

Whitening with \bm{z} = \mathbf{S}^{-1}(\bm{x} - \bm{\mu}) turns this into a quadratic in standard normal \bm{z} with

\begin{aligned} \widetilde{\mathbf{Q}}_2 &= \mathbf{S}\mathbf{Q}_2\mathbf{S}, \\ \tilde{\bm{q}}_1 &= 2\mathbf{S}\mathbf{Q}_2\bm{\mu} + \mathbf{S}\bm{q}_1, \\ \tilde{q}_0 &= q(\bm{\mu}). \end{aligned} \qquad [\text{eq. 2}]

1. Quadratic domains: generalized chi-square §2.1

1. If \widetilde{\mathbf{Q}}_2 = 0, the boundary is flat and the probability is \Phi(\tilde{q}_0 / \lVert \tilde{\bm{q}}_1 \rVert), with \Phi the standard normal cdf §2.1.
2. Otherwise, eigen-decompose \widetilde{\mathbf{Q}}_2 = \mathbf{R}\mathbf{D} \mathbf{R}', with \mathbf{R} orthogonal, and set \bm{b} = \mathbf{R}' \tilde{\bm{q}}_1. The quadratic becomes a sum of scaled non-central chi-squares with one degree of freedom each, plus a normal term §2.1.
3. Merge terms with equal eigenvalue w_j (multiplicity k_j) and map to the generalized chi-square parameters \bm{w}, \bm{k}, \bm{\lambda}, s, m §2.1:

\lambda_j = \frac{1}{4 w_j^2} \sum_{i:\, D_i = w_j} b_i^2, \quad s = \sqrt{\sum_{i'} b_{i'}^2}, \quad m = q(\bm{\mu}) - \bm{w}\cdot\bm{\lambda}

Here i indexes nonzero eigenvalues and i' the zero ones §2.1.

4. The required probability is then a one-dimensional integral of the generalized chi-square cdf at 0. The paper names Ruben's, Imhof's, IFFT and ray methods for it §2.1.

2. Arbitrary domains: ray-trace §2.2

Standardise \tilde{f}(\bm{z}) = f(\mathbf{S}\bm{z} + \bm{\mu}) §2.2.1. Write each standard point as z\bm{n}, with \bm{n} a unit direction and z the coordinate along it. The integral becomes an integral over directions of an axial integral §2.2.1.

Along one ray, the domain is described by two things §2.2.2:

• the initial sign \psi(\bm{n}) at z = -\infty: 1 inside, -1 outside, 0 grazing;
• the crossing points z_i(\bm{n}), the zeros of the ray function.

For a quadratic domain with coefficients \tilde{q}_2(\bm{n}) = \bm{n}' \widetilde{\mathbf{Q}}_2\bm{n}, \tilde{q}_1(\bm{n}) = \tilde{\bm{q}}_1'\bm{n} and \tilde{q}_0, the three cases are combined as §2.2.2:

\psi(\bm{n}) = \operatorname{sign}\!\left[ 4\,\operatorname{sign}(\tilde{q}_2) - 2\,\operatorname{sign}(\tilde{q}_1) + \operatorname{sign}(\tilde{q}_0) \right] \quad [\text{§2.2.2}]

The standard normal along a ray has the complementary distribution §2.2.3

\bar{\Phi}^{\text{ray}}_k(z) = \frac{1 - \operatorname{sign}(z)\, F_{\chi^2_k}(z^2)}{2},

where F_{\chi^2_k} is the chi-squared cdf with k degrees of freedom. The total solid angle is \Omega_k, equal to 2 in 1D, 2\pi in 2D, 4\pi in 3D and 2\pi^2 in 4D §2.2.3.

The probability in the angular slice d\bm{n} is §2.2.4

dp(\bm{n}) = \frac{\alpha(\bm{n})}{\Omega_k}\, d\bm{n}, \qquad \alpha(\bm{n}) = \psi(\bm{n}) + 1 + 2\psi(\bm{n}) \sum_i (-1)^i\, \bar{\Phi}^{\text{ray}}_k\!\left(z_i(\bm{n})\right).

The sum runs over the roots along the ray, in increasing order §2.2.4.

# Probability of one ray slice  [§2.2.4]
# Input:  initial_sign in {1, 0, -1}; roots sorted ascending
# Output: alpha (scalar)

def ray_alpha(initial_sign, roots, dim):
    # Start from the sign at -inf, shifted to {0, 1, 2}  [§2.2.4]
    alpha = initial_sign + 1
    # Each crossing adds or removes the tail beyond it  [§2.2.4]
    for index, root in enumerate(roots, start=1):
        alpha += 2 * initial_sign * (-1) ** index * ray_survival(root, dim)
    return alpha


# Survival of the standard normal along a ray  [§2.2.3]
def ray_survival(z, dim):
    return (1 - sign(z) * chi_squared_cdf(z ** 2, dim)) / 2

Integration over directions §2.2.4:

• Up to 4D the paper integrates \alpha on a grid spanning half the angular space, using any standard quadrature. An adaptive grid can follow sharp changes of the boundary §2.2.4.
• Beyond 4D, or in general, it draws random directions. It samples the standard multinormal and normalises the vectors to get uniform directions \bm{n}. The estimate is \langle \alpha(\bm{n}) \rangle / 2 §2.2.4.

# Ray-trace probability over a set of directions  [§2.2.4]
# Input:  directions (m, dim) unit vectors covering half the sphere
#         weights (m,): quadrature weights on that half-sphere
# Output: p (scalar)

def ray_trace_probability(directions, weights, trace, dim):
    total = 0
    for direction, weight in zip(directions, weights):
        # Sign and crossings along this direction  [§2.2.2]
        initial_sign, roots = trace(direction)
        total += weight * ray_alpha(initial_sign, roots, dim)
    # Divide by the total solid angle  [§2.2.3] [§2.2.4]
    return total / omega(dim)

For Monte Carlo, the same \alpha is averaged over uniform directions on the full sphere, and the estimate is that mean divided by 2 §2.2.4. Set operations (complement, intersection, union) on domains are applied to the ray traces directly §2.2.5. For a function f(\bm{x}) of a normal vector, the cdf is the normal probability of c - f(\bm{x}) > 0, and the pdf is obtained by numerical differentiation of that cdf §2.2.6.

3. Classification between two normals §3.1

Classes a and b have priors p_a, p_b and value matrix \mathbf{V} with entries v_{ij} §3. The Bayes decision variable \beta(\bm{x}) for two normals is the quadratic eq. 3:

\begin{aligned} \beta(\bm{x}) &= \bm{x}'\mathbf{Q}_2\bm{x} + \bm{q}_1'\bm{x} + q_0, \\ \mathbf{Q}_2 &= \tfrac{1}{2}\left(\bm{\Sigma}_b^{-1} - \bm{\Sigma}_a^{-1}\right), \\ \bm{q}_1 &= \bm{\Sigma}_a^{-1}\bm{\mu}_a - \bm{\Sigma}_b^{-1}\bm{\mu}_b, \\ q_0 &= \tfrac{1}{2}\left(\bm{\mu}_b'\bm{\Sigma}_b^{-1}\bm{\mu}_b - \bm{\mu}_a'\bm{\Sigma}_a^{-1}\bm{\mu}_a + \ln\frac{|\bm{\Sigma}_b|}{|\bm{\Sigma}_a|}\right) + \ln\frac{p_a v_a}{p_b v_b}. \end{aligned} \qquad [\text{eq. 3}]

Here v_a = v_{aa} - v_{ab} and v_b = v_{bb} - v_{ba} §3 (the paper defines v_i = v_{ii} - \sum_{j\ne i} v_{ij} for the general case §3). Deciding a when \beta > 0 maximises the expected value gain §3.1. With \mathbf{V} = \mathbf{1} it minimises the overall error §3.1.

Single interval (yes/no) §3.1.1. For 1D normals with \sigma_a > \sigma_b, the Bayes error rates are eq. 4:

\begin{aligned} p(B \mid a) &= P\!\left(\chi'^2_{1,\sigma_a^2\lambda} < \sigma_b^2 c\right), \qquad p(A \mid b) = P\!\left(\chi'^2_{1,\sigma_b^2\lambda} > \sigma_a^2 c\right), \\ \lambda &= \left(\frac{\mu_a - \mu_b}{\sigma_a^2 - \sigma_b^2}\right)^2, \qquad c = \lambda + \frac{2\ln(\sigma_a/\sigma_b)}{\sigma_a^2 - \sigma_b^2}. \end{aligned}

The overall error p(e) is the average of the two §3.1.1.

Two-interval task §3.1.2. Two stimuli come from a and b. The optimal rule compares the two log likelihood ratios l_{a/b}(x_1) and l_{a/b}(x_2), not the stimuli. The pair log ratio is l_{ab/ba} = l_{a/b}(x_1) - l_{a/b}(x_2) §3.1.2. For normals, the boundary is two perpendicular lines when the variances differ §3.1.2 fig. 3c.

m-interval task §3.1.3. One of m stimuli is the signal. The maximum accuracy is

p(c) = \int_{-\infty}^{\infty} F_b^{\,m-1}(l)\, f_a(l)\, dl,

where f_a and F_b are the pdf and cdf of the log likelihood ratio l §3.1.3. The text uses Wickens' eq. 6.19 for this result §3.1.3.

Discriminability indices §3.1.4. For Bayes error p(e) with equal priors, the Bayes index is

d'_b = -2\,Z\!\left(p(e)\right) = 2\,Z\!\left(p(c)\right),

where Z is the standard normal quantile §3.1.4. The other indices are §3.1.4:

• Simpson and Fitter: d'_a = \lVert \mathbf{S}_{\text{rms}}^{-1} \bm{d}\rVert with \mathbf{S}_{\text{rms}} = [(\bm{\Sigma}_a + \bm{\Sigma}_b)/2]^{1/2};
• Egan and Clarke: d'_e = \lVert \mathbf{S}_{\text{avg}}^{-1} \bm{d}\rVert with \mathbf{S}_{\text{avg}} = (\mathbf{S}_a + \mathbf{S}_b)/2;
• here \bm{d} = \bm{\mu}_a - \bm{\mu}_b.

For two bivariate normals with correlation \rho eq. 5:

d'^2 = \frac{d_x'^2 + d_y'^2 - 2\rho\, d'_x d'_y}{1 - \rho^2}, \qquad d'_x = \frac{\mu_{b,x} - \mu_{a,x}}{\sigma_x}, \quad d'_y = \frac{\mu_{b,y} - \mu_{a,y}}{\sigma_y}.

The sign of the mean difference is kept in d'_x and d'_y §3.1.4.

ROC curves §3.1.5. Sweeping a criterion on the likelihood ratio l traces the ROC curve. d'_b is the accuracy at the point that maximises hit minus false-alarm rate §3.1.5.

Custom boundaries §3.1.6. The optimal boundary can be replaced by a user-supplied quadratic, a ray trace or an implicit function §3.1.6.

4. Classification from samples §3.2 to §3.6

• Fit §3.2: maximum-likelihood means, covariances and priors are the sample means, sample covariances and relative counts. The optimal boundary follows from these §3.2.
• Sample-optimal boundary \gamma §3.2: start from \beta and optimise its (k+1)(k+2)/2 free quadratic parameters to maximise the outcome value of the given samples. The paper says this stays in the quadratic family and limits overfitting §3.2.

# Sample-optimal quadratic boundary  [§3.2]
# Input:  samples_a (n_a, k), samples_b (n_b, k), beta (quadratic)
# Output: gamma (quadratic with (k+1)(k+2)/2 free parameters)

def sample_optimal_boundary(samples_a, samples_b, beta):
    # Flatten the upper triangle of Q2, then q1 and q0  [§3.2]
    params = flatten_quadratic(beta)
    # Maximise the outcome value of the samples, starting from beta  [§3.2]
    params = maximise(sample_outcome_value, params, samples_a, samples_b)
    return unflatten_quadratic(params)

• Scaling discriminability §3.3. Move the signal distribution b towards the noise distribution a by a factor r, interpolating the mean and the sd matrix (not the covariance):

# Pull signal b towards noise a by linear interpolation  [§3.3]
# r = 1 leaves the pair unchanged; r -> 0 makes b match a
# The sd matrices sd_a, sd_b are symmetric square roots of covariances

def squeeze_signal(mu_a, sd_a, mu_b, sd_b, r):
    mu_b_new = r * mu_b + (1 - r) * mu_a
    sd_b_new = r * sd_b + (1 - r) * sd_a
    return mu_b_new, sd_b_new

For samples, each point of b is whitened, then mapped by the new parameters §3.3. The paper also scales the decision variable l towards its median l_m by l' = l - l_m(1 - r) §3.3.

• Multiple normals §3.4. Each class region is the intersection of its pairwise quadratic regions, f(\bm{x}) = \min_i q_n^i(\bm{x}) > 0. This is not quadratic, so only the ray method computes its error rates §3.4.
• Dimension reduction §3.5. A decision variable on a group of dimensions replaces that group. Projecting onto one dimension makes the integral or the classification visible in 1D §3.5.
• Normality check §3.6. Compare observed and predicted hit and false-alarm counts along a family of boundaries, such as changing the prior or scaling the covariances. Counts are binomial, or multinomial for several classes. The paper uses 100 trials as typical of a psychophysics experiment §3.6 fig. 7b.

Hyperparameters

The paper states few numerical settings. Values below come from the paper's text or figures, or from code defaults. Code defaults are not in the paper.

Results

Author-reported. Setup: MATLAB implementations on an AMD Ryzen Threadripper 2950X 16-core machine, timings to 1% precision, Monte Carlo as the reference §5 fig. 8b. In the table, "0" means the difference lies within the spread of the Monte Carlo runs §5.

Source for all rows: fig. 8b. The table sits inside the benchmark figure float in the source.

The paper reports the normal-model checks and applications qualitatively in addition to these numbers §6.

Paper vs code

Code: MIT License, code_commit 4929451. The repository is by the first author. Checked by reading the source only; no code was run.

Agrees with the paper:

• Standardisation of the quadratic standard_quad.m matches eq. 2.
• The initial-sign formula standard_ray_trace.m matches the combined \psi of §2.2.2.
• Roots of a quadratic along a ray, with tangents removed roots_ray_quad.m, match "crossing points only" §2.2.2.
• The symbolic ray probability prob_ray_sym.m matches \alpha(\bm{n}) §2.2.4: the chi-squared survival (1 - \operatorname{sign}(z) F)/2 and the alternating sum.
• The two-normal quadratic opt_class_quad.m matches eq. 3.
• d'_e and d'_a classify_normals.m match §3.1.4.
• Discriminability scaling, normal and sample paths [code: classify_normals.m#L294-L300, L316-L323], matches §3.3.

Differs from the paper:

1. Grid versus Monte Carlo cutoff. The paper uses a grid for angles up to 4D §2.2.4. The code defaults to grid integration up to 3D and to Monte Carlo above, and force_mc overrides the grid integrate_normal.m classify_normals.m.
2. Quadratic method switch. For quadratic domains above 3D the default is gx2 integrate_normal.m. If gx2 returns 0 and the method was defaulted, the ray method is used instead integrate_normal.m. The paper states no cutoff.
3. Sample-optimal boundary. The paper starts from \beta and optimises the outcome value §3.2. In optimize_samp_bd.m the default is 'svm', which fits a degree-2 polynomial SVM and returns that boundary without starting from \beta optimize_samp_bd.m optimize_samp_bd.m. The doc warns the SVM boundary can do worse than \beta classify_normals.m. The paper's approach corresponds to the fminsearch and multi-start branches optimize_samp_bd.m. classify_normals sets samp_opt to 100 by default classify_normals.m, but it does not pass that value to the optimiser classify_normals.m, so the SVM is used unless the caller sets it.
4. d'_b outside equal priors. The paper defines d'_b = -2Z(p(e)) with equal priors §3.1.4. The code uses that formula only when priors are 0.5 and no custom domain is given classify_normals.m classify_normals.m. Otherwise it uses Z(\text{hit}) - Z(\text{false alarm}) from the error matrix classify_normals.m. The two agree only when the two error rates are equal. The doc describes this value as "Bayes-optimal" classify_normals.m, although with a custom boundary it is the value for that boundary.
5. Tiny errors. When p(e) underflows, the code uses an asymptotic d'_b \approx 2\sqrt{w - \ln w} classify_normals.m. The paper describes no asymptotic for multi-dimensional cases. It says d'_b becomes infinite beyond about 75 and that d'_e can replace it in 1D only §3.1.4 §5. The code falls back to d'_e in 1D classify_normals.m.
6. Implicit-domain roots. The paper says any algorithm that finds all zeros on an interval can be used §2.2.2. The code instead samples fun_resol (default 100) equally spaced points on [-\text{fun\_span}, +\text{fun\_span}], looks for sign changes and refines each with fzero ray_trace_fun.m. See Pitfalls.
7. Missing integrator in this snapshot. int_norm_ray and Phibar_ray_split are called but not defined in the repository at this commit prob_ray.m integrate_normal.m. The angular integration, the tolerance logic, the GPU batching and the split tail sums therefore cannot be checked here. The GPU speed-up claim §2.2.4 is also unverified.
8. gx2 functions are external. norm_quad_to_gx2_params and gx2cdf are not in this repository; they come from the separate gx2 toolbox the paper names int_norm_quad_gx2.m §2.1. Not fetched, so not verified.

Not stated in the paper: whether sample covariances use n or n-1. The code uses MATLAB cov, which divides by n-1 classify_normals.m.

Pitfalls

• Finite search radius. A ray method for non-quadratic domains only sees crossings within radius m. The error is the normal mass beyond radius m, which in k dimensions is the chi tail \bar{F}_{\chi_k}(m) §2.2.2. The paper gives no numeric m. The code uses 5, and the figures trace to 3 sd integrate_normal.m fig. 2a.
• Starting sign. The paper defines \psi at z = -\infty §2.2.2. The code reads the sign at -\text{fun\_span} ray_trace_fun.m. The two agree only if the domain does not change sign beyond the span.
• Missed roots. The sign-change scan misses two crossings inside one grid step and any tangency. With the default resolution the step is 10/99 \approx 0.10 sd ray_trace_fun.m integrate_normal.m. Tangency is handled for quadratics only roots_ray_quad.m. Increase fun_resol or check the domain by hand for thin features.
• Tolerance. The default RelTol is 10^{-2} integrate_normal.m. Answers are 1% accurate unless the caller tightens it.
• Monte Carlo for tiny probabilities. The benchmark shows sampling fails beyond d' \approx 10 §5 fig. 8a. Ray and quadratic methods stay accurate to about d' = 75, the double-precision limit §5. Beyond that, p(e) = 0 and d'_b = \infty in the paper's double-precision form §5. The code's log10 encoding and asymptote extend this differently (see Paper vs code, item 5).
• Rounding. The paper says it keeps tiny and large summands separate to avoid rounding and uses accurate tail probabilities §5. The code's prob_ray.m falls back to a symbolic evaluation only when vpa=true and the numeric result is zero prob_ray.m. The numeric split itself is in the missing helper (see Paper vs code, item 7).
• Sensitivity of d'_a. d'_a \le d'_e \le d'_b in 1D. d'_a can understate discriminability by about 30% and is not a good Bayes index §3.1.4.
• d'_b is not a metric. It is symmetric and positive but breaks the triangle inequality §3.1.4.
• Two-interval rule. The optimal rule is not "pick the larger stimulus" unless the variances are equal §3.1.2.
• Interpolate sd matrices, not covariances. Points are linear in the mean and the sd matrix, so only those should be blended §3.3.
• Normality. The check in §3.6 compares predicted and observed counts for a given number of trials. The paper leaves it to the investigator to judge whether a deviation is small enough for their problem §3.6 fig. 7b.
• Dependencies. The gx2 path needs an external toolbox that is not in this repository int_norm_quad_gx2.m.

Lineage

• Builds on: Ruben's probability-content results §2.1; Imhof's quadratic-form method §2.1; Genz and Bretz's ray-style integration for rectangles §2.2; Genz's ellipsoid code §2.1; Simpson and Fitter's d'_a and Egan and Clarke's d'_e §3.1.4; Wickens' m-interval accuracy §3.1.3. None of these has a card here, so they are named only.
• Later work by the same first author: Das (2024), arXiv:2404.05062, which the paper says gives the inverse map from generalized chi-square parameters to quadratic forms, analytic pdfs, and log or variable-precision arithmetic for d' beyond 75 §2.1 §2.2.6 §5.
• Compared against: Monte Carlo sampling §5; d'_a and d'_e §3.1.4.
• Superseded by: not known from content already in this wiki.