Monte Carlo Integration

Monte Carlo integration uses random sampling to estimate integrals — particularly powerful in high dimensions where deterministic quadrature becomes infeasible.

Basic Idea

To compute $\int_a^b f(x)\, dx$, draw $N$ uniform random samples $x_i \in [a, b]$ and estimate:

$\int_a^b f(x)\, dx \approx \frac{b-a}{N} \sum_{i=1}^N f(x_i)$

The error scales as $1/\sqrt{N}$ regardless of dimension — a huge advantage over grid-based methods in $d$ dimensions (which scale as $1/N^{2/d}$).

Importance Sampling

When $f$ is sharply peaked, plain sampling is inefficient. Instead, sample from a distribution $p(x)$ that is large where $f$ is large:

$\int f(x)\, dx = \int \frac{f(x)}{p(x)}\, p(x)\, dx \approx \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{p(x_i)}$

Optimal choice: $p(x) \propto |f(x)|$ (reduces variance to zero for positive $f$).

Estimating π

A classic application: generate random points $(x, y) \in [0,1)^2$. The fraction landing inside the unit quarter-circle ($x^2 + y^2 < 1$) converges to $\pi/4$:

$\pi \approx \frac{4 \cdot \text{(points inside circle)}}{N}$

Reproducible Randomness: LCG

For reproducibility, we use a Linear Congruential Generator (LCG):

$s_{n+1} = (a \cdot s_n + c) \bmod m, \qquad x_n = s_n / m$

with parameters $a = 1{,}664{,}525$, $c = 1{,}013{,}904{,}223$, $m = 2^{32}$ (Numerical Recipes constants). Starting from a seed $s_0$, this generates a deterministic sequence of uniform pseudo-random numbers in $[0, 1)$.

Implementation

• lcg_random(seed, n) — returns list of $n$ uniform $[0,1)$ values via LCG
• monte_carlo_integrate(f_func, a, b, N=10000, seed=42) — plain Monte Carlo estimate
• monte_carlo_pi(N=10000, seed=42) — estimate $\pi$ using the unit circle method

import math

def lcg_random(seed, n):
    a = 1664525
    c = 1013904223
    m = 2**32
    state = seed
    result = []
    for _ in range(n):
        state = (a * state + c) % m
        result.append(state / m)
    return result