3. One-Dimensional Monte Carlo Integration

Ray Tracing: The Rest of Your Life (v3.2.3): 3 One Dimensional MC Integration / 3.3 One-Dimensional Monte Carlo Integration

Section 3.2 derived an estimate of $\pi$ from an area ratio. I now generalize the same idea into a framework for estimating arbitrary integrals with the Monte Carlo method. This framework forms the foundation for the rest of Book 3.

This chapter covers two topics:

• The general form of a Monte Carlo estimator using uniform sampling
• Variance reduction through importance sampling with a carefully chosen PDF

The General Form of Monte Carlo Integration

Using a random variable $x\sim p(x)$, an integral can be expressed as an expected value:

$$I=\int f(x)dx=\int \frac{f(x)}{p(x)}p(x)dx={\mathbb{E}}_{x\sim p}\left[\frac{f(x)}{p(x)}\right]$$

A Monte Carlo estimator approximates this expectation with a sample mean:

$${\hat{I}}_N=\frac{1}{N}\sum _{i=1}^N\frac{f(x_i)}{p(x_i)},x_i\sim p(x)$$

For any distribution $p(x)$, the expected value of ${\hat{I}}_N$ equals $I$, making it an unbiased estimator. The choice of $p$ affects the variance.

Estimation with Uniform Sampling

Consider estimating the following integral with the uniform distribution $p(x)=1/2$ on $[0,2]$:

$$I={\int }_0^2{x}^{2}dx$$

The exact value is $I=8/3$. The estimator becomes

$${\hat{I}}_N=\frac{1}{N}\sum _{i=1}^N\frac{x_i^2}{1/2}=\frac{2}{N}\sum _{i=1}^N{x}_i^2$$

Importance Sampling

Reducing variance requires reducing the variation in $f(x)/p(x)$. Drawing more samples where $f(x)=x^2$ is large—that is, at larger values of $x$—makes the information contributed by each sample more uniform.

Choose $p(x)=x/2$ on $[0,2]$. This distribution follows the same increasing trend as $f(x)=x^2$. Its cumulative distribution function (CDF) is

$$P(x)={\int }_0^x\frac{t}{2}dt=\frac{x^2}{4}$$

Using inverse transform sampling, a sample can be drawn from $u\sim U[0,1)$ as

$$x=P^{-1}(u)=\sqrt{4u}$$

The estimator is

$${\hat{I}}_N=\frac{1}{N}\sum _{i=1}^N\frac{x_i^2}{x_i/2}$$

Ideally, $p(x)\propto f(x)$, but even a distribution that merely follows the same trend can have lower variance than uniform sampling.

Rust Implementation

The r303-one-dimensional-mc crate is another standalone implementation with no dependency on common. It contains the same local XorShift64 implementation as Section 3.2.

The PDF and inverse-transform sampling functions are as follows.

r303-one-dimensional-mc/src/lib.rs

fn uniform_pdf(_x: f64) -> f64 {
    0.5
}

fn linear_pdf(x: f64) -> f64 {
    0.5 * x
}

fn sample_uniform(rng: &mut XorShift64) -> f64 {
    rng.range_f64(0.0, 2.0)
}

fn sample_linear_pdf(rng: &mut XorShift64) -> f64 {
    let u = rng.next_f64();
    (4.0 * u).sqrt()
}

Each estimation function accumulates $f(x)/p(x)$ over $N$ samples and returns the mean.

r303-one-dimensional-mc/src/lib.rs

fn estimate_integral_uniform(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let x = sample_uniform(rng);
        sum += (x * x) / uniform_pdf(x);
    }
    sum / samples as f64
}

fn estimate_integral_importance(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let x = sample_linear_pdf(rng);
        sum += (x * x) / linear_pdf(x);
    }
    sum / samples as f64
}

Because a single run is sensitive to random chance, summarize_trials repeats the same conditions and calculates the mean absolute error. Separate RNG instances for the uniform and importance-sampling estimators prevent the trials for one method from affecting the other.

r303-one-dimensional-mc/src/lib.rs

fn summarize_trials(samples: usize, trials: usize, seed: u64) -> (f64, f64) {
    let mut rng_u = XorShift64::new(seed ^ 0x3030_a1u64);
    let mut rng_i = XorShift64::new(seed ^ 0x3030_b2u64);

    let mut sum_abs_error_uniform = 0.0;
    let mut sum_abs_error_importance = 0.0;

    for _ in 0..trials {
        let estimate_u = estimate_integral_uniform(samples, &mut rng_u);
        let estimate_i = estimate_integral_importance(samples, &mut rng_i);
        sum_abs_error_uniform += (estimate_u - EXACT_INTEGRAL).abs();
        sum_abs_error_importance += (estimate_i - EXACT_INTEGRAL).abs();
    }

    (
        sum_abs_error_uniform / trials as f64,
        sum_abs_error_importance / trials as f64,
    )
}

generate_report constructs its output in two stages. It first compares the two methods in a single estimate with $N={10}^6$, then records the mean absolute error across repeated trials in CSV format for smaller sample counts from $N=100$ through $10000$.

Differences Between C++ and Rust

The original C++ uses the global random-number state provided by rand().

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

Rust passes rng explicitly as an argument. Because summarize_trials maintains rng_u and rng_i as independent instances, the trial sequences for the two methods do not interfere with each other and remain reproducible.

r303-one-dimensional-mc/src/lib.rs

const EXACT_INTEGRAL: f64 = 8.0 / 3.0;

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 uniform_pdf(_x: f64) -> f64 {
    0.5
}

fn linear_pdf(x: f64) -> f64 {
    0.5 * x
}

fn sample_uniform(rng: &mut XorShift64) -> f64 {
    rng.range_f64(0.0, 2.0)
}

fn sample_linear_pdf(rng: &mut XorShift64) -> f64 {
    let u = rng.next_f64();
    (4.0 * u).sqrt()
}

fn estimate_integral_uniform(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let x = sample_uniform(rng);
        sum += (x * x) / uniform_pdf(x);
    }
    sum / samples as f64
}

fn estimate_integral_importance(samples: usize, rng: &mut XorShift64) -> f64 {
    let mut sum = 0.0;
    for _ in 0..samples {
        let x = sample_linear_pdf(rng);
        sum += (x * x) / linear_pdf(x);
    }
    sum / samples as f64
}

fn summarize_trials(samples: usize, trials: usize, seed: u64) -> (f64, f64) {
    let mut rng_u = XorShift64::new(seed ^ 0x3030_a1u64);
    let mut rng_i = XorShift64::new(seed ^ 0x3030_b2u64);

    let mut sum_abs_error_uniform = 0.0;
    let mut sum_abs_error_importance = 0.0;

    for _ in 0..trials {
        let estimate_u = estimate_integral_uniform(samples, &mut rng_u);
        let estimate_i = estimate_integral_importance(samples, &mut rng_i);
        sum_abs_error_uniform += (estimate_u - EXACT_INTEGRAL).abs();
        sum_abs_error_importance += (estimate_i - EXACT_INTEGRAL).abs();
    }

    (
        sum_abs_error_uniform / trials as f64,
        sum_abs_error_importance / trials as f64,
    )
}

pub fn generate_report() -> String {
    let mut output = String::new();
    output.push_str("One Dimensional Monte Carlo Integration\n");
    output.push_str("=====================================\n\n");

    output.push_str("Target integral: I = integral_0^2 x^2 dx = 8/3\n");
    output.push_str(&format!("Exact value: {:.12}\n\n", EXACT_INTEGRAL));

    let samples = 1_000_000usize;
    let mut rng_uniform = XorShift64::new(0x3030_0001u64);
    let mut rng_importance = XorShift64::new(0x3030_0002u64);

    let uniform_est = estimate_integral_uniform(samples, &mut rng_uniform);
    let importance_est = estimate_integral_importance(samples, &mut rng_importance);

    output.push_str("Single run (N=1,000,000):\n");
    output.push_str(&format!(
        "method,estimate,abs_error\nuniform,{:.12},{:.12}\nimportance_x_over_2,{:.12},{:.12}\n\n",
        uniform_est,
        (uniform_est - EXACT_INTEGRAL).abs(),
        importance_est,
        (importance_est - EXACT_INTEGRAL).abs()
    ));

    output.push_str("Average absolute error across repeated trials:\n");
    output.push_str("samples,trials,uniform,importance_x_over_2\n");

    let settings = [(100usize, 200usize), (1_000usize, 200usize), (10_000usize, 100usize)];
    for (trial_samples, trials) in settings {
        let (err_u, err_i) = summarize_trials(trial_samples, trials, 0x3030_9000u64 + trial_samples as u64);
        output.push_str(&format!(
            "{},{},{:.12},{:.12}\n",
            trial_samples, trials, err_u, err_i
        ));
    }

    output
}

Results in the Browser

The method,estimate,abs_error rows compare the two methods in a single estimate with $N={10}^6$. In each samples,trials,uniform,importance_x_over_2 row, the first two columns give the trial conditions and the last two give the mean absolute errors. The difference is more pronounced at smaller sample counts, where importance_x_over_2 is usually lower.

Summary

• I derived the general Monte Carlo estimator ${\hat{I}}_N=\frac{1}{N}\sum f(x_i)/p(x_i)$. It is unbiased for any $p$, while the choice of $p$ affects its variance.
• I estimated ${\int }_0^2{x}^{2}dx$ using uniform sampling with $p(x)=1/2$.
• Importance sampling with $p(x)=x/2$ can be generated by the inverse transform $x=\sqrt{4u}$ and reduces variance for the same sample count.
• By using summarize_trials to calculate mean absolute error over multiple trials, I measured the difference in variance that can be difficult to observe in a single run.

The next chapter extends this $f(x)/p(x)$ structure to integration over a sphere.