2. A Simple Monte Carlo Program

Ray Tracing: The Rest of Your Life (v3.2.3): 2 A Simple Monte Carlo Program / 3.2 A Simple Monte Carlo Program

The defining feature of the Monte Carlo method (MC method) is that a simple program can statistically approximate a true value, with accuracy improving as it runs longer. This chapter uses the estimation of $\pi$ as a minimal example and establishes the theoretical foundation for all of Book 3.

This chapter covers three topics:

• Constructing a Monte Carlo estimator from random samples
• Empirically confirming the $O(1/\sqrt{N})$ convergence rate
• Reducing variance through stratified sampling

Estimating π

Generate uniformly random points $(x,y)$ in the square $[-1,1]\times [-1,1]$ and count the proportion that fall inside the unit circle $x^2+y^2<1$. The square has area $4$ and the unit circle has area $\pi$, so

$$\frac{\text{points inside}}{N}\approx \frac{\pi }{4}\Rightarrow \pi \approx 4\frac{\text{points inside}}{N}$$

This directly exhibits the structure of replacing the expected value of a random variable with a sample mean, making it the simplest example of Monte Carlo integration.

The estimate approaches the true value as the sample count $N$ increases, but the convergence rate of $O(1/\sqrt{N})$ is slow. Multiplying $N$ by $100$ improves accuracy by only a factor of $10$.

Rust Implementation

Because this chapter prepares the probability theory independently of ray-tracing geometry, the r302-monte-carlo crate is a standalone implementation with no dependency on common.

To ensure reproducibility with a fixed seed, I implement a local Xorshift64 random number generator. An equivalent implementation exists in common::utils, but because the internals of common have not been introduced at this point, the full implementation appears here.

r302-monte-carlo/src/lib.rs

struct XorShift64 {
    state: u64,
}

impl XorShift64 {
    fn new(seed: u64) -> Self {
        let state = if seed == 0 { 0x9e3779b97f4a7c15 } else { seed };
        Self { state }
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
    }

    fn range_f64(&mut self, min: f64, max: f64) -> f64 {
        min + (max - min) * self.next_f64()
    }
}

The estimation function draws $N$ samples and estimates $\pi$ from the number of points inside the circle.

r302-monte-carlo/src/lib.rs

fn estimate_pi(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut inside_circle = 0usize;
    for _ in 0..samples {
        let x = rng.range_f64(-1.0, 1.0);
        let y = rng.range_f64(-1.0, 1.0);
        if x * x + y * y < 1.0 {
            inside_circle += 1;
        }
    }
    4.0 * inside_circle as f64 / samples as f64
}

generate_report constructs its output in two stages. It first uses a single estimate with N=1000 to show the accuracy of one trial, then records convergence in CSV format at checkpoints from $N=10$ through ${10}^6$. The third block, added in the next section, compares stratified sampling.

r302-monte-carlo/src/lib.rs

pub fn generate_report() -> String {
    let mut output = String::new();
    output.push_str("Monte Carlo Estimation of PI\n");
    output.push_str("============================\n\n");

    let mut first_run_rng = XorShift64::new(0x30200001);
    let first_estimate = estimate_pi(1000, &mut first_run_rng);
    output.push_str(&format!(
        "Single run (N=1000): {:.12} (error={:.12})\n\n",
        first_estimate,
        (first_estimate - PI).abs()
    ));

    output.push_str("Convergence checkpoints:\n");
    output.push_str("samples,estimate,abs_error\n");

    let checkpoints = [10usize, 100, 1_000, 10_000, 100_000, 1_000_000];
    let mut checkpoint_index = 0usize;
    let mut inside_circle = 0usize;
    let mut rng = XorShift64::new(0x30200002);

    for n in 1..=checkpoints[checkpoints.len() - 1] {
        let x = rng.range_f64(-1.0, 1.0);
        let y = rng.range_f64(-1.0, 1.0);
        if x * x + y * y < 1.0 {
            inside_circle += 1;
        }

        if n == checkpoints[checkpoint_index] {
            let estimate = 4.0 * inside_circle as f64 / n as f64;
            output.push_str(&format!(
                "{n},{estimate:.12},{error:.12}\n",
                n = n,
                estimate = estimate,
                error = (estimate - PI).abs()
            ));
            checkpoint_index += 1;
            if checkpoint_index == checkpoints.len() {
                break;
            }
        }
    }

    // ...stratified sampling comparison (next section)...
    output
}

Stratified Sampling

The slow $O(1/\sqrt{N})$ convergence arises from the accidental clustering inherent in purely random samples. Stratified sampling eliminates this clustering by dividing the sample space into equal cells and drawing exactly one sample from each cell.

Divide $[-1,1{]}^2$ into a $\sqrt{N}\times \sqrt{N}$ grid and choose one point from cell $(i,j)$ using uniform random values $(u,v)\in [0,1{)}^2$:

$$x=2\frac{i+u}{\sqrt{N}}-1,y=2\frac{j+v}{\sqrt{N}}-1$$

Because every cell contains exactly one point, the samples cover the domain more uniformly and their variance decreases. Unstratified samples can accidentally cluster in some regions and leave others empty; stratification prevents this.

The final block of generate_report compares unstratified and stratified sampling with the same total sample count: $\sqrt{N}=1000$, or ${10}^6$ samples. Each $(i,j)$ iteration draws one unstratified sample and one stratified sample, keeping the comparison conditions equal.

r302-monte-carlo/src/lib.rs

    let sqrt_n: usize = 1000;
    let total = sqrt_n * sqrt_n;
    let mut inside_regular = 0usize;
    let mut inside_stratified = 0usize;
    let mut strat_rng = XorShift64::new(0x30200003);

    for i in 0..sqrt_n {
        for j in 0..sqrt_n {
            // Regular: unconstrained uniform sample in [-1, 1)^2
            let xr = strat_rng.range_f64(-1.0, 1.0);
            let yr = strat_rng.range_f64(-1.0, 1.0);
            if xr * xr + yr * yr < 1.0 {
                inside_regular += 1;
            }
            // Stratified: one jittered sample per (i, j) cell
            let xs = 2.0 * ((i as f64 + strat_rng.next_f64()) / sqrt_n as f64) - 1.0;
            let ys = 2.0 * ((j as f64 + strat_rng.next_f64()) / sqrt_n as f64) - 1.0;
            if xs * xs + ys * ys < 1.0 {
                inside_stratified += 1;
            }
        }
    }

    let est_regular = 4.0 * inside_regular as f64 / total as f64;
    let est_stratified = 4.0 * inside_stratified as f64 / total as f64;
    let err_regular = (est_regular - PI).abs();
    let err_stratified = (est_stratified - PI).abs();
    output.push_str(&format!("\nStratified sampling comparison (sqrt_N={sqrt_n}):\n"));
    output.push_str(&format!(
        "regular    estimate={est_regular:.12} error={err_regular:.12}\n"
    ));
    output.push_str(&format!(
        "stratified estimate={est_stratified:.12} error={err_stratified:.12}\n"
    ));

Even with the same ${10}^6$ samples, stratified sampling tends to produce a smaller error. This advantage diminishes as dimensionality increases, however, because of the curse of dimensionality. A $\sqrt{N}\times \sqrt{N}$ partition in two dimensions becomes $N^{1/3}\times N^{1/3}\times N^{1/3}$ in three, sharply reducing how many subdivisions each dimension receives for the same total sample count. Practical ray tracers commonly apply stratification only to low-dimensional components such as pixel sampling and use ordinary random numbers inside high-dimensional integrals.

Differences Between C++ and Rust

The original C++ generates random numbers from the global state used by rand().

inline double random_double() {
    return rand() / (RAND_MAX + 1.0);
}

Rust has no global random-number state equivalent to rand(). Because this chapter does not depend on common, it implements Xorshift64 locally and passes rng explicitly as an argument. Restricting a function's stateful effects to its arguments and return value also improves reproducibility.

For output formatting, C++ uses std::cout << std::fixed << std::setprecision(12), whereas Rust specifies the format directly inside the string with format!("{:.12}").

r302-monte-carlo/src/lib.rs

const PI: f64 = std::f64::consts::PI;

struct XorShift64 {
    state: u64,
}

impl XorShift64 {
    fn new(seed: u64) -> Self {
        let state = if seed == 0 { 0x9e3779b97f4a7c15 } else { seed };
        Self { state }
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
    }

    fn range_f64(&mut self, min: f64, max: f64) -> f64 {
        min + (max - min) * self.next_f64()
    }
}

fn estimate_pi(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut inside_circle = 0usize;
    for _ in 0..samples {
        let x = rng.range_f64(-1.0, 1.0);
        let y = rng.range_f64(-1.0, 1.0);
        if x * x + y * y < 1.0 {
            inside_circle += 1;
        }
    }
    4.0 * inside_circle as f64 / samples as f64
}

pub fn generate_report() -> String {
    let mut output = String::new();
    output.push_str("Monte Carlo Estimation of PI\n");
    output.push_str("============================\n\n");

    // Single run
    let mut first_run_rng = XorShift64::new(0x30200001);
    let first_estimate = estimate_pi(1000, &mut first_run_rng);
    output.push_str(&format!(
        "Single run (N=1000): {:.12} (error={:.12})\n\n",
        first_estimate,
        (first_estimate - PI).abs()
    ));

    // Convergence checkpoints
    output.push_str("Convergence checkpoints:\n");
    output.push_str("samples,estimate,abs_error\n");

    let checkpoints = [10usize, 100, 1_000, 10_000, 100_000, 1_000_000];
    let mut checkpoint_index = 0usize;
    let mut inside_circle = 0usize;
    let mut rng = XorShift64::new(0x30200002);

    for n in 1..=checkpoints[checkpoints.len() - 1] {
        let x = rng.range_f64(-1.0, 1.0);
        let y = rng.range_f64(-1.0, 1.0);
        if x * x + y * y < 1.0 {
            inside_circle += 1;
        }

        if n == checkpoints[checkpoint_index] {
            let estimate = 4.0 * inside_circle as f64 / n as f64;
            output.push_str(&format!(
                "{n},{estimate:.12},{error:.12}\n",
                n = n,
                estimate = estimate,
                error = (estimate - PI).abs()
            ));
            checkpoint_index += 1;
            if checkpoint_index == checkpoints.len() {
                break;
            }
        }
    }

    // Stratified sampling comparison
    let sqrt_n: usize = 1000;
    let total = sqrt_n * sqrt_n;
    let mut inside_regular = 0usize;
    let mut inside_stratified = 0usize;
    let mut strat_rng = XorShift64::new(0x30200003);

    for i in 0..sqrt_n {
        for j in 0..sqrt_n {
            // Regular: unconstrained uniform sample in [-1, 1)^2
            let xr = strat_rng.range_f64(-1.0, 1.0);
            let yr = strat_rng.range_f64(-1.0, 1.0);
            if xr * xr + yr * yr < 1.0 {
                inside_regular += 1;
            }
            // Stratified: one jittered sample per (i, j) cell
            let xs = 2.0 * ((i as f64 + strat_rng.next_f64()) / sqrt_n as f64) - 1.0;
            let ys = 2.0 * ((j as f64 + strat_rng.next_f64()) / sqrt_n as f64) - 1.0;
            if xs * xs + ys * ys < 1.0 {
                inside_stratified += 1;
            }
        }
    }

    let est_regular = 4.0 * inside_regular as f64 / total as f64;
    let est_stratified = 4.0 * inside_stratified as f64 / total as f64;
    let err_regular = (est_regular - PI).abs();
    let err_stratified = (est_stratified - PI).abs();
    output.push_str(&format!("\nStratified sampling comparison (sqrt_N={sqrt_n}):\n"));
    output.push_str(&format!(
        "regular    estimate={est_regular:.12} error={err_regular:.12}\n"
    ));
    output.push_str(&format!(
        "stratified estimate={est_stratified:.12} error={err_stratified:.12}\n"
    ));

    output
}

Results in the Browser

At the convergence checkpoints, abs_error does not decrease monotonically. It fluctuates with random clustering but becomes smaller overall. The sampling comparison shows that stratification tends to produce a smaller error with the same ${10}^6$ samples, although a particular random sequence can reverse the result.

Summary

• I implemented a Monte Carlo estimator using the minimal example of estimating $\pi$. Its convergence rate is $O(1/\sqrt{N})$, so multiplying $N$ by $100$ improves accuracy by only a factor of $10$.
• Convergence checkpoints make the relationship between sample count and accuracy observable.
• Stratified sampling covers the domain more uniformly with the same number of samples and can reduce variance in two dimensions. The benefit diminishes in higher dimensions because of the curse of dimensionality.

The next chapter takes this estimate of $\pi$ as its starting point and develops the general framework of Monte Carlo integration, treating the integral of $f(x)$ as the mean of a random variable.