Quasi-Monte Carlo Integration

Why Use Quasi-Monte Carlo?

Likelihoods for models with latent variables often require integration over quantities that are not observed. When the integral has no convenient analytic form, simulation-based integration provides a numerical approximation. Ordinary Monte Carlo uses pseudorandom draws; quasi-Monte Carlo (QMC) instead uses points chosen to cover the integration region evenly. This regular coverage can reduce integration error, particularly when the integrand is reasonably smooth and the effective dimension is modest.

The method is not guaranteed to outperform pseudorandom simulation for every integrand or sample size. Its value depends on the shape of the integrand, the way latent variables are transformed, and the dimension of the integration problem. A useful implementation therefore treats QMC as a numerical design whose accuracy must be assessed, rather than as a guarantee of precision.

From a Likelihood Integral to a Unit Cube

Suppose an observed outcome \(y\) has conditional likelihood \(p(y\mid z,\theta)\) given a latent variable \(Z=z\), whose density is \(f(z\mid\theta)\). The likelihood contribution after integrating out the latent variable is

If \(Z\) has a continuous cumulative distribution function \(F_\theta\), the probability integral transform gives \(Z=F_\theta^{-1}(U)\) for \(U\) uniform on \([0,1]\). Thus the same contribution can be written as an integral over a fixed unit interval. For a \(d\)-dimensional latent vector, a suitable transformation \(T_\theta\) similarly gives

The transformation must reproduce the latent-variable distribution in the original integral. If the transformation instead defines a proposal distribution, the corresponding likelihood weights must also be included; omitting them changes the integral. The earlier material on importance sampling for likelihood integrals treats that weighted alternative.

Discrepancy and the QMC Error Principle

For a point set, star discrepancy measures how far its empirical coverage of axis-aligned boxes \([0,t_1]\times\cdots\times[0,t_d]\) differs from the volume of those boxes. Low discrepancy means that this difference is small over the collection of such boxes. It captures a geometric property of the points, not a measure of the likelihood integrand itself.

The Koksma–Hlawka inequality connects that geometric property to integration error. For a function \(h\) of bounded Hardy–Krause variation \(V_{\mathrm{HK}}(h)\), it states

where \(D_N^*\) is the star discrepancy of the point set. The inequality explains the basic QMC idea: an evenly covering design can give a small error when the transformed integrand also has manageable variation. For fixed dimension, well-known low-discrepancy sequences have discrepancy of order \((\log N)^d/N\), up to constants. This is a theoretical bound, not usually a practical error estimate: the variation can be hard to calculate, and the bound can be loose. In higher dimensions, the constants and the factor involving dimension can matter greatly.

Worked Likelihood Integrals

Randomized QMC and Error Assessment

A deterministic QMC estimate is reproducible, but it does not provide a standard Monte Carlo standard error from the variation of independent draws: the points were not independent draws. A common solution is randomized QMC. A valid randomization, such as a suitable scramble of a low-discrepancy sequence, preserves its space-filling structure while making the randomized estimate unbiased. Generate \(R\) independent randomizations, each with \(N\) points, and compute one estimate \(\widehat I_r\) from each.

If each randomized point is uniform on the unit cube, linearity of expectation gives the unbiasedness justification:

The points within one randomized design can be dependent; independence between the \(R\) randomizations is what supports estimating the variation across replicate estimates. Let \(\bar I_R\) be the mean of those estimates and \(s_R\) their sample standard deviation. The estimated standard error of their mean is \(s_R/\sqrt{R}\). This measures numerical variability under the chosen randomized design. It does not measure uncertainty from sampling the observed data.

How QMC Compares With Pseudorandom Integration

With pseudorandom independent uniform draws \(U_j\), the usual estimate is also an average, but its standard error is estimated from the sample variation of \(h(U_j)\), divided by \(\sqrt{N}\). Under finite variance, its error typically decreases at rate \(N^{-1/2}\). QMC instead aims to improve coverage of the integration region; for sufficiently regular integrands and fixed dimension, its error can decrease faster as \(N\) grows. This does not make QMC universally more accurate at a fixed \(N\). A discontinuous or sharply concentrated integrand, an awkward transformation, or a high effective dimension can weaken its advantage.

When comparing methods, use the same integral, transformation, and computational budget where possible. Randomized QMC replicates allow a numerical variability estimate, while independent pseudorandom replications provide a comparable assessment for Monte Carlo. A single unusually favorable estimate is not enough to establish that one method is more accurate. It is also useful to inspect whether increasing \(N\) stabilizes the likelihood contribution and whether estimates from independent randomizations agree.

Common Pitfall and Practical Checks

Before relying on a QMC likelihood integral, verify that the transformation represents the intended latent distribution and that any required weights are present. Then compare estimates across independent randomizations and, when practical, across increasing point counts. If the estimates remain unstable, QMC has not yet demonstrated adequate numerical accuracy; more points, a better transformation, or a different integration strategy may be needed.