Alternatives

The probabilities this package computes are multivariate normal rectangle and argmax probabilities, and good software for those predates this project by decades. This page lists the implementations people already use, says when each is the right choice, and describes the reference improvements shipped as winning.alternatives — which are new, and therefore not as battle-tested as the incumbents they modify. Every comparison number below comes from a committed script under research/alternatives/.

Existing implementations, by language

R

mvtnorm is the stalwart: pmvnorm computes one rectangle probability by Genz–Bretz randomized quasi–Monte Carlo, documented for dimensions up to 1000 (up to 20 for the exact Miwa algorithm), at roughly three hundred thousand downloads a month. Recent versions add lpRR and slpRR: reduced-rank rectangle probabilities and their scores for covariance of the form BB′ + diagonal, the Marsaglia–Genz–Bretz construction. For one specified winner under factor covariance, lpRR computes essentially the same integral this engine does — the difference appears when every winner is wanted at once, which is the case this package exists for. TruncatedNormal implements Botev's minimax exponential tilting, the strongest single rare-orthant estimator we know. tlrmvnmvt combines Genz separation-of-variables with tile-low-rank covariance algebra for very high-dimensional single rectangles. mlogit estimates multinomial probit by GHK behind its familiar interface. VeccTMVN reaches linear cost for a single region by combining Vecchia dependence truncation with minimax tilting, and a 2026 doubly-blocked Genz implementation reports order-of-magnitude wall-time gains over the packages above for large single rectangles through cache blocking, SIMD and threading. The pattern across all of them: one region per call, so pricing all N winner events costs N calls however fast each call is.

Python

scipy.stats.multivariate_normal.cdf wraps Genz's quasi–Monte Carlo routine for rectangle probabilities. For the full multinomial-probit choice vector there is no standard Python tool; the scattered research implementations of GHK on GitHub are unmaintained. This package fills that slot, and winning.alternatives adds the strongest GHK-family construction we could build for cross-checking it.

Stata

GHK's native home: Mata ships ghk() and ghkfast() (the latter caches its simulation points across likelihood iterations), driving cmmprobit, asmprobit and the community mvprobit. One reading note: the factor(#) option parameterizes the covariance as low-rank-plus-diagonal, but the evaluator still Cholesky-factorizes each difference covariance and runs dense GHK — the factor structure changes the model, not the computation.

Our suggested improvements

Both ship in winning.alternatives, both are new, and neither has the decades of production hardening the packages above carry. Certify anything important against Monte Carlo argmax frequencies; the module docstrings say how.

The reduced-rank bridge

reduced_rank_representation(mu, V, D, i) returns the Marsaglia rectangle representation of one winner's probability — the (k+1)-column loading matrix, residual diagonal, and bounds — ready for lpRR or any rectangle solver, and verified against pmvnorm in this repository's tests. The complete per-winner protocol also ships compiled (per_winner_reduced_rank_shares, a Rust kernel with a numpy fallback), so the wall-clock comparison below is same-toolchain. Use it when you want a second opinion on a single share from machinery independent of this package. Pricing the complete vector this way costs one call per winner: measured at n = 1000 (rank 2, common scrambled-Sobol draws, both sides compiled), eighteen times the engine's exact 0.125 seconds at total-variation 0.036, or 148 times at 0.0045, with the per-winner cost growing as the square of the field size and deep-tail winners falling below the estimator's floor at finite sample counts.

The CDF-gradient evaluator

cdf_gradient_shares(mu, V, D) computes every winner probability from one rectangle estimator: the identity is that the share vector is the integrated boundary gradient of the maximum's CDF along the diagonal, the estimator is GHK with the factors as a Kalman state (no dense Cholesky), and one reverse-mode sweep per grid point returns all n components, so the whole vector is linear in n. We found no published GHK protocol combining these pieces; as far as we know this is the strongest GHK-family all-share evaluator in existence, and we give it away because the shared field still dominates it: conditional on the factors, its draws simulate idiosyncratic dimensions the field integrates analytically. Measured at n = 200 against a 131k-node reference: at ranks 4, 8 and 16 the engine's quasi–Monte Carlo factor nodes are five to sixty-five times more accurate at equal or less time, and at rank 2 the margin is fifty-fold and up. The evaluator remains useful as an independent cross-check, a GPU starting point, and an honest baseline. It requires jax, and one warning transfers: its output sums to one identically for any estimate of the CDF, so total mass is not a diagnostic of correctness.

Which to use

One rectangle in high dimension with general covariance: mvtnorm, or TruncatedNormal when the probability is rare, or tlrmvnmvt when the dimension is in the thousands. One winner's share under factor covariance: lpRR, or the bridge above. The complete share vector, its Jacobian, removal counterfactuals, or calibration back from observed shares, under factor, block or tree covariance: that combination is what winning.factor is for, and the convergence demo shows the comparison live.