Monte Carlo Integration

Prerequisites

• Linear Algebra — matrix operations (for understanding high-dimensional sampling)
• Basic probability theory — random variables, expectation, variance
• Basic calculus — the definition of the definite integral and change of variables

The following probability-theory background is recommended:

Basic Concepts

From Definite Integrals to Expectations

Consider the definite integral

$$I = \int_\Omega f(\mathbf{x}), d\mathbf{x}$$

where $\Omega \subseteq \mathbb{R}^d$ is the integration domain. The core idea of Monte Carlo methods is to rewrite the integral as an expectation.

Let $\mathbf{X}$ be a random variable uniformly distributed over $\Omega$, with probability density function

$$p(\mathbf{x}) = \frac{1}{V}, \qquad V = \text{vol}(\Omega)$$

Then the integral becomes

$$I = V \cdot \mathbb{E}[f(\mathbf{X})]$$

The Basic Estimator

The basic Monte Carlo estimator: draw $N$ independent samples $\mathbf{x}_1, \dots, \mathbf{x}_N$ from $p(\mathbf{x})$ and define

$$\hat{I} = \frac{V}{N} \sum_{i=1}^{N} f(\mathbf{x}_i)$$

Properties:

1. Unbiasedness: $\mathbb{E}[\hat{I}] = I$ (its expectation equals the true value exactly).
2. Variance: $\text{Var}(\hat{I}) = \frac{V^2}{N} \text{Var}(f(\mathbf{X}))$.
3. Convergence rate: the standard error of $\hat{I}$ is $\sigma_{\hat{I}} \propto 1/\sqrt{N}$, independent of the dimension $d$.

The last property is the core advantage of Monte Carlo methods: the error of classical numerical integration (trapezoidal rule, Simpson's rule) grows exponentially with the dimension (the "curse of dimensionality"), whereas the convergence rate of Monte Carlo is always $O(N^{-1/2})$.

The Law of Large Numbers and the Central Limit Theorem

Law of large numbers (Kolmogorov's strong law). If $X_1, X_2, \dots$ are independent identically distributed random variables with $\mathbb{E}[|X_1|] < \infty$, then

$$\bar{X}n = \frac{1}{n}\sum{i=1}^{n} X_i ;\xrightarrow{\text{a.s.}}; \mathbb{E}[X_1]$$

Applied to Monte Carlo: $\hat{I}/V \to \mathbb{E}[f(\mathbf{X})]$ holds almost surely, hence $\hat{I} \to I$.

Central limit theorem (CLT). Under the above conditions, if $\text{Var}(X_1) = \sigma^2 < \infty$, then

$$\frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}} ;\xrightarrow{d}; \mathcal{N}(0, 1)$$

Applied to Monte Carlo: for sufficiently large $N$, $\hat{I}$ is approximately normally distributed:

$$\hat{I} ;\sim; \mathcal{N}!\left(I,; \frac{V^2 \sigma_f^2}{N}\right)$$

where $\sigma_f^2 = \text{Var}(f(\mathbf{X}))$. This provides the basis for confidence intervals.

Variance Estimation and Confidence Intervals

In practice $\sigma_f^2$ is unknown and is estimated by the sample variance:

$$\hat{\sigma}f^2 = \frac{1}{N-1} \sum{i=1}^{N} \bigl(f(\mathbf{x}i) - \bar{f}\bigr)^2, \qquad \bar{f} = \frac{1}{N}\sum{i=1}^{N} f(\mathbf{x}_i)$$

The standard error of $\hat{I}$ is estimated as

$$\hat{\sigma}_{\hat{I}} = \frac{V}{\sqrt{N}} \hat{\sigma}_f$$

The $95%$ confidence interval is $\hat{I} \pm 1.96, \hat{\sigma}_{\hat{I}}$.

Generalization to weighted sampling: when the sampling density $p(\mathbf{x})$ is not uniform (importance sampling), each sample $i$ carries the weight $w_i = 1/p(\mathbf{x}_i)$, and

$$\hat{I} = \frac{1}{N} \sum_{i=1}^{N} w_i \cdot f(\mathbf{x}_i)$$

the variance estimate becomes

$$\hat{\sigma}^2 = \frac{1}{N} \sum_{i=1}^{N} w_i^2 f(\mathbf{x}i)^2 - \left(\frac{1}{N}\sum{i=1}^{N} w_i f(\mathbf{x}_i)\right)^2$$

oCAS's StatisticsAccumulator uses a similar weighted-variance form (see the "Statistics Accumulator" section below).

Core Theory

Importance Sampling

The variance of a Monte Carlo estimator depends on how much the integrand "fluctuates" under the sampling distribution. Importance sampling reduces the variance by choosing a non-uniform sampling distribution $p(\mathbf{x})$.

Principle: rewrite the integral as

$$I = \int_\Omega \frac{f(\mathbf{x})}{p(\mathbf{x})} \cdot p(\mathbf{x}), d\mathbf{x} = \mathbb{E}_p!\left[\frac{f(\mathbf{X})}{p(\mathbf{X})}\right]$$

The estimator is

$$\hat{I}p = \frac{1}{N} \sum{i=1}^{N} \frac{f(\mathbf{x}_i)}{p(\mathbf{x}_i)}, \qquad \mathbf{x}_i \sim p$$

Optimal sampling distribution: the variance-minimizing $p^*$ is

$$p^*(\mathbf{x}) = \frac{|f(\mathbf{x})|}{\int_\Omega |f(\mathbf{y})|, d\mathbf{y}}$$

when $f$ does not change sign, $w \cdot f = \int |f| \cdot \mathrm{sgn}(f)$ is constant and the variance is zero; if $f$ changes sign the variance is still minimized. But $p^$ itself requires knowing the integral — a chicken-and-egg problem. In practice, an adaptive grid is used to approximate $p^$.

Key constraint: $p(\mathbf{x})$ must be strictly positive at every point where $f(\mathbf{x}) \neq 0$. If $p$ vanishes where $f$ is non-zero, the estimator is unbounded and the variance is infinite.

The Vegas Adaptive Grid

The Vegas algorithm (Lepage, 1978) is the most widely used adaptive Monte Carlo integration method. Its core idea is to approximate the optimal distribution $p^*$ by a piecewise-constant sampling density.

The Product-Grid Structure

Vegas works on the $d$-dimensional unit hypercube $[0,1]^d$ (the physical domain is mapped to the hypercube by a linear change of variables). The sampling density has a product form:

$$p(\mathbf{x}) = \prod_{k=1}^{d} p_k(x_k)$$

Each $p_k$ is an independent one-dimensional piecewise-constant density on the $k$-th axis. This means the grid of each dimension is maintained independently — the $d$-dimensional problem is decomposed into $d$ one-dimensional problems.

The grid of the $k$-th dimension: divide $[0,1]$ into $M$ bins with boundaries

$$0 = b_0 < b_1 < \cdots < b_M = 1$$

The width of the $j$-th bin is $\Delta_j = b_{j+1} - b_j$. On that bin $p_k(x) = 1/(M \cdot \Delta_j)$ (constant), because each bin is chosen with equal probability.

The Sampling Procedure

For each sample point $\mathbf{x} = (x_1, \dots, x_d)$:

1. Pick a bin: for the $k$-th dimension, choose the bin number $j_k$ uniformly at random from ${0, 1, \dots, M-1}$.
2. Pick a point: choose $x_k$ uniformly at random within $[b_{j_k}, b_{j_k+1}]$.
3. Compute the Jacobian (weight): the probability density at this point is $p_k(x_k) = 1/(M \cdot \Delta_{j_k})$, hence

$$w_k = \frac{1}{p_k(x_k)} = M \cdot \Delta_{j_k}$$

The total weight is

$$w = \prod_{k=1}^{d} w_k = \prod_{k=1}^{d} (M \cdot \Delta_{j_k})$$

4. Accumulate: add $w \cdot f(\mathbf{x})$ to the statistics accumulator.

Physical meaning of the weight: $w = 1/p(\mathbf{x})$ is the reciprocal of the sampling density — samples in sparse regions (large $\Delta_j$) get larger weights and samples in dense regions (small $\Delta_j$) get smaller weights, compensating for the non-uniform sampling.

Grid Adaptation: The Training Signal

In each iteration, Vegas records the "importance" of every bin. For the $j$-th bin of the $k$-th dimension, the training signal is

$$d_j = \sum_{\text{samples in bin } j} w \cdot f^2(\mathbf{x})$$

Intuition: $f^2$ measures the "energy" (variance contribution) of the integrand in the region of the bin, and $w$ compensates for non-uniform sampling. Bins with high $d_j$ contribute heavily to the integration error and need a finer grid.

Grid Updating: Smoothed Cumulative-Arc-Length Redistribution

Given the training signals ${d_0, d_1, \dots, d_{M-1}}$, Vegas updates the grid boundaries through the following steps:

Step 1: 3-bin smoothing. For each bin $j$, compute the smoothed value:

$$\tilde{d}j = \frac{d{j-1} + d_j + d_{j+1}}{3}$$

Missing neighbours of boundary bins ($j = 0$ or $j = M-1$) are replaced by 0. The purpose of smoothing is to suppress the noise of individual bins and make the grid changes more stable.

Step 2: mean normalization. Compute the average $\bar{d} = \frac{1}{M}\sum_j \tilde{d}_j$ and set

$$\hat{d}_j = \frac{\tilde{d}_j}{\bar{d}}$$

After normalization, $\hat{d}_j = 1$ means "average importance" and $\hat{d}_j > 1$ means the bin needs a finer grid.

Step 3: damping. Introduce a learning rate $\alpha > 0$ (default $\alpha = 1.5$ in oCAS) and apply a power transform to the normalized values:

$$\hat{d}_j ;\leftarrow; \hat{d}_j^{1/\alpha}$$

When $\alpha = 1$ there is no damping ($\hat{d}_j^{1/1} = \hat{d}_j$); when $\alpha > 1$ the exponent is $< 1$, compressing extreme values and making grid updates more conservative. This prevents the grid from overreacting to the noise of a single iteration.

Step 4: cumulative arc length and equidistant redistribution. Compute the cumulative sums:

$$S_j = \sum_{i=0}^{j-1} \hat{d}_i, \qquad j = 0, 1, \dots, M$$

where $S_0 = 0$ and $S_M = \sum_{i=0}^{M-1} \hat{d}_i$. Treat $S$ as an "arc-length" curve; the new bin boundaries are equidistributed along the arc length:

$$b'_j = \text{the position in } [0,1] \text{ corresponding to the arc length } S = j \cdot S_M / M$$

Concretely, for $j = 1, \dots, M-1$, find the index $i$ in $S$ with $S_i \leq j \cdot S_M / M < S_{i+1}$ and then interpolate linearly:

$$b'j = \frac{i + \frac{j \cdot S_M / M - S_i}{S{i+1} - S_i}}{M}$$

Fix $b'_0 = 0$ and $b'_M = 1$.

Step 5: monotonicity clamping. Due to floating-point precision and numerical noise, the $b'_j$ may not be strictly increasing. Repair the monotonicity of the result:

$$b'_j ;\leftarrow; \max(b'j,; b'{j-1})$$

Intuition: the geometric meaning of arc-length redistribution is that — if we plot the cumulative-importance curve in the $(x, S(x))$ plane, the "equal-arc-length division points" make every new bin sweep equal arc length on the curve, i.e. every new bin carries approximately equal integration contribution.

Stratified Sampling

Stratified sampling is another variance-reduction technique. Its idea is to divide the integration domain into disjoint subdomains ("strata") and sample independently within each stratum:

$$I = \sum_{j=1}^{M} I_j, \qquad I_j = \int_{\Omega_j} f(\mathbf{x}), d\mathbf{x}$$

The estimator of each stratum is $\hat{I}_j = V_j \cdot \bar{f}_j$ ($V_j$ the stratum volume, $\bar{f}_j$ the within-stratum sample mean), and the total estimator is $\hat{I} = \sum_j \hat{I}_j$.

Variance analysis:

$$\text{Var}(\hat{I}) = \sum_{j=1}^{M} \frac{V_j^2}{N_j} \sigma_j^2$$

where $N_j$ is the number of samples in the $j$-th stratum and $\sigma_j^2$ is the variance of $f$ within the $j$-th stratum. If the stratum variances differ, the optimal allocation (Neyman allocation) is

$$N_j \propto V_j \sigma_j$$

i.e. strata with higher variance receive more samples.

The relation between Vegas and stratified sampling: the product grid of Vegas is essentially an adaptive stratified sampling — the bin boundaries stratify the integration domain and sampling is uniform within each stratum. The difference is that Vegas's stratum boundaries are adapted iteratively rather than fixed in advance.

Inverse-Variance Weighting Across Iterations

Vegas runs several iterations (with $N$ samples each); every iteration produces an estimate $\hat{I}_r$ and its standard error $\hat{\sigma}_r$. How should the iteration results be combined?

Inverse-variance weighting: let $w_r = 1/\hat{\sigma}_r^2$; the combined estimate is

$$\hat{I}_{\text{combined}} = \frac{\sum_r w_r \hat{I}_r}{\sum_r w_r}$$

and the combined error is

$$\hat{\sigma}_{\text{combined}} = \frac{1}{\sqrt{\sum_r w_r}}$$

Optimality: when the estimates of different iterations are mutually independent, inverse-variance weighting is the minimum-variance unbiased way of combining them (a direct corollary of the Gauss–Markov theorem).

Online inverse-variance accumulation: to avoid storing the estimates of all iterations, StatisticsAccumulator uses online updates (a weighted mean plus an incremental $\chi^2$; this is not Welford's algorithm itself). It maintains the state variables:

• integral: the current combined estimate
• error: the current combined standard error
• chi_square: the across-iteration $\chi^2$ statistic

When iteration $r$ completes (estimate $\mu_r$, error $\sigma_r$), the update rules are:

$$w_{\text{prev}} = 1/\sigma_{\text{prev}}^2, \quad w_r = 1/\sigma_r^2, \quad w_{\text{new}} = w_{\text{prev}} + w_r$$

$$\hat{I}{\text{new}} = \frac{w{\text{prev}} \cdot \hat{I}{\text{prev}} + w_r \cdot \mu_r}{w{\text{new}}}$$

$$\sigma_{\text{new}} = 1/\sqrt{w_{\text{new}}}$$

Implementation detail: to prevent zero-variance iterations from overflowing the weights, oCAS clamps the error at $\sigma_r \geq 10^{-150}$ ($10^{-300}$ is still representable as an f64).

The $\chi^2$ Diagnostic

Purpose: to test whether the estimates of the different iterations are consistent — if the $\hat{I}_r$ differ far more than their error bars, the variance is underestimated or the integrator has not converged.

Definition:

$$\chi^2 = \sum_{r=1}^{R} w_r (\hat{I}r - \hat{I}{\text{combined}})^2$$

where $w_r = 1/\hat{\sigma}_r^2$.

Interpretation: if the iterations are independent and the variance estimates are accurate, $\chi^2$ follows a chi-square distribution with $R - 1$ degrees of freedom. The empirical rule is:

oCAS's StatisticsAccumulator updates $\chi^2$ online at the end of each iteration:

$$\chi^2_{\text{new}} = \chi^2_{\text{prev}} + w_{\text{prev}} (\hat{I}{\text{prev}} - \hat{I}{\text{new}})^2 + w_r (\mu_r - \hat{I}_{\text{new}})^2$$

Note this is an incremental update — the estimates of past iterations need not be stored.

Implementation in oCAS

oCAS's Vegas implementation lives in ocas-eval/src/numeric/vegas.rs; the accompanying statistics accumulator is in ocas-eval/src/numeric/statistics.rs.

Data Structures

Vegas
├── opts: VegasOptions        // tunable options
├── axes: Vec<GridAxis>       // one independent grid per dimension
│   ├── boundaries: Vec<f64>  // M+1 bin boundaries ∈ [0,1]
│   └── bin_accum: Vec<f64>   // training-signal accumulator for the M bins
└── accumulator: StatisticsAccumulator  // inverse-variance weighting across iterations
    ├── sum_w, sum_wf, sum_wf2  // weighted sums for the current iteration
    ├── integral, error          // combined estimate and error
    ├── chi_square               // χ² diagnostic
    └── iterations               // number of completed iterations

Default Parameters

Implementing the Sampling Procedure

Pseudocode of one iteration:

#![allow(unused)]
fn main() {
for _ in 0..n_samples {
    // 1. sample each dimension independently
    let mut x = Vec::with_capacity(n_dims);
    let mut jac = 1.0;
    for axis in axes.iter_mut() {
        // pick a bin uniformly
        let b = rng.random_range(0..n_bins);
        let lo = boundaries[b];
        let hi = boundaries[b + 1];
        // sample uniformly within the bin
        let u = rng.random::<f64>();
        let xi = lo + (hi - lo) * u;
        // Jacobian = M × bin width = 1/pdf
        let wi = (hi - lo) * n_bins as f64;
        x.push(xi);
        jac *= wi;
    }
    // 2. evaluate and accumulate
    let fx = f(&x);
    accumulator.add_sample(jac, fx);
    // 3. record the training signal
    for (i, xi) in x.iter().enumerate() {
        axes[i].add_training(*xi, jac, fx * fx);
    }
}
}

Key implementation details:

• RNG choice: uses Xoshiro256PlusPlus (rather than thread_rng) with a fixed seed to guarantee reproducible results.
• Product Jacobian: the total weight is jac = ∏ w_k — the dimensions are sampled independently, so the total density is the product of the per-dimension densities.
• Training signal: add_training locates the bin by binary search ($O(\log M)$) and accumulates weight × f² into the corresponding bin.

Implementing the Grid Update

The implementation steps of GridAxis::update(learning_rate):

#![allow(unused)]
fn main() {
fn update(&mut self, learning_rate: f64) {
    let n = self.bin_accum.len();
    let total: f64 = self.bin_accum.iter().sum();
    if total <= 0.0 { return; }

    // 1. 3-bin smoothing
    let avg = total / n as f64;
    let mut d = vec![0.0; n];
    for i in 0..n {
        let prev = if i > 0 { bin_accum[i-1] } else { 0.0 };
        let next = if i+1 < n { bin_accum[i+1] } else { 0.0 };
        d[i] = (prev + bin_accum[i] + next) / 3.0 / avg;
    }

    // 2. damping: d[i] = d[i]^(1/learning_rate)
    if learning_rate != 1.0 {
        for v in d.iter_mut() {
            *v = v.max(1e-30).powf(1.0 / learning_rate);
        }
    }

    // 3. cumulative arc length
    let mut cum = vec![0.0; n+1];
    for i in 0..n { cum[i+1] = cum[i] + d[i]; }

    // 4. equidistant redistribution
    let mut new_boundaries = vec![0.0; n+1];
    new_boundaries[0] = 0.0;
    new_boundaries[n] = 1.0;
    let mut j = 0;
    for i in 1..n {
        let target = i as f64 / n as f64 * cum[n];
        while j < n && cum[j+1] < target { j += 1; }
        let frac = (target - cum[j]) / (cum[j+1] - cum[j]);
        new_boundaries[i] = (j as f64 + frac) / n as f64;
    }

    // 5. monotonicity clamping
    for i in 1..=n {
        if new_boundaries[i] < new_boundaries[i-1] {
            new_boundaries[i] = new_boundaries[i-1];
        }
    }
    new_boundaries[n] = 1.0;
    self.boundaries = new_boundaries;
    self.bin_accum.fill(0.0);  // reset the training signal
}
}

Numerical safeguards in the implementation:

• the normalised values take max(1e-30) before the power transform, guaranteeing $d_i > 0$ so the cumulative arc length is strictly increasing (preventing division by zero in the interpolation).
• when learning_rate = 1.0, the power transform (together with the 1e-30 clamp) is skipped, so cum[j+1] - cum[j] can still be zero (a bin completely flat), in which case frac = 0.
• the final clamping guarantees the boundaries are strictly non-decreasing, with the endpoints fixed at 0 and 1.

The Statistics Accumulator

StatisticsAccumulator accumulates online within each iteration using three variables:

$$S_w = \sum_i w_i, \qquad S_{wf} = \sum_i w_i f_i, \qquad S_{wf^2} = \sum_i w_i f_i^2$$

At the end of an iteration, the estimates are

$$\mu = \frac{S_{wf}}{S_w}, \qquad \sigma^2 = \frac{S_{wf^2}}{S_w} - \mu^2$$

Note this is not the standard sample-variance formula — it is the weighted-variance formula, applicable when every sample has a different weight (the Jacobian weights of Vegas).

The One-Dimensional Convenience Function

integrate_1d(f, a, b, opts) internally performs the linear change of variables $x = a + u \cdot (b-a)$:

#![allow(unused)]
fn main() {
let width = b - a;
let wrapped = |u: &[f64]| f(a + u[0] * width) * width;
let mut vegas = Vegas::new(1, opts);
vegas.integrate(&wrapped)
}

The Jacobian factor $(b-a)$ is multiplied directly into the integrand; the user-supplied $f(x)$ receives physical coordinates rather than hypercube coordinates.

Advanced Topics

When to Use Vegas vs. Symbolic Integration

Limitations of Vegas

1. Product-grid assumption: Vegas assumes the optimal sampling density factorizes as $\prod_k p_k(x_k)$. When the integrand has strong correlations between variables (such as $f(x,y) = \delta(x - y)$), the product grid cannot approximate it effectively and convergence slows down.

2. Piecewise-constant approximation: the grid density is piecewise constant — sampling is uniform within each bin. If the integrand varies sharply inside a bin (e.g. a narrow peak), the variance estimate for that bin may be inaccurate.

3. Adaptation bias: early iterations use a grid that has not yet converged, so their samples are not fully independent. Vegas mitigates this via inverse-variance weighting — converged iterations get larger weights.

4. Deterministic seed: oCAS uses a fixed RNG seed (default 0x0C45) for reproducibility. This is useful in tests and benchmarks, but production environments should consider a random seed to avoid systematic bias.

Possible Improvements

• Vegas+ (Lepage, 2020): introduces a "reflection" strategy into the grid update to handle symmetric integrands.
• Suave (Hahn, 2005): combines the advantages of the Vegas adaptive grid with stratified sampling.
• Divonne (Friedman & Harris, 1996): uses partitioning trees and heuristic subdivision strategies.
• Cuhre: purely deterministic quadrature (not Monte Carlo), suitable for low-dimensional smooth integrands.

These methods could be added to oCAS by implementing the Integrator trait.

References

1. Lepage, G. P. "A new algorithm for adaptive multidimensional integration." Journal of Computational Physics, 27(2):192–203, 1978. — The original paper on the Vegas algorithm.

2. James, F. "Monte Carlo theory and practice." Reports on Progress in Physics, 43(9):1145–1189, 1980. — A survey of Monte Carlo methods, including a detailed description of Vegas.

3. Lepage, G. P. "Vegas: An adaptive multi-dimensional integration program." Cornell preprint CLNS 80-447, 1980. — Implementation details and usage guide for Vegas.

4. Hahn, T. "CUBA—a library for multidimensional numerical integration." Computer Physics Communications, 168(2):78–95, 2005. — A comparison of Vegas, Suave, Divonne, and Cuhre.

5. Neal, R. M. "Annealed importance sampling." Statistics and Computing, 11(2):125–139, 2001. — The theoretical foundations of importance sampling.

6. Kahn, H. & Marshall, A. W. "Methods of reducing sample size in Monte Carlo computations." Journal of the Operations Research Society of America, 1(5):263–278, 1953. — Early work on stratified sampling and importance sampling.