3. 一次元モンテカルロ積分

Ray Tracing: The Rest of Your Life (v3.2.3): 3 One Dimensional MC Integration / 3.3 一次元モンテカルロ積分

3.2 では $\pi$ の推定を「面積比」として導きました。ここでは同じアイデアを一般化し，任意の積分をモンテカルロ法で推定する枠組みを作ります。この枠組みが第 3 編全体を貫く基礎となります。

扱うのは 2 つのトピックです。

• 一様サンプリングによるモンテカルロ推定量の一般形
• PDF を工夫した重点サンプリング（importance sampling）による分散低減

モンテカルロ積分の一般形

確率変数 $x\sim p(x)$ を使うと，積分を期待値として書けます。

$$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]$$

この期待値を標本平均で近似するのがモンテカルロ推定量です。

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

$p(x)$ がどんな分布でも ${\hat{I}}_N$ の期待値は $I$ に等しく，不偏推定量です。$p$ の選び方は分散に影響します。

一様サンプリングによる推定

$[0,2]$ 上の一様分布 $p(x)=1/2$ を使って

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

を推定します。厳密値は $I=8/3$ です。推定量は

$${\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$$

となります。

重点サンプリング

分散を下げるには，$f(x)/p(x)$ の変動を小さくすればよいです。$f(x)=x^2$ が大きい領域（$x$ が大きいところ）でより多くサンプルを取れば，1 サンプルあたりの情報量が均一になります。

$[0,2]$ で $p(x)=x/2$ を選びます。これは $f(x)=x^2$ に比例した分布です。累積分布（CDF）は

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

逆関数法により，$u\sim U[0,1)$ から

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

でサンプルできます。推定量は

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

となります。$p(x)\propto f(x)$ が理想ですが，完全に比例しなくても傾向が同じであれば一様よりも分散が小さくなります。

Rust 実装

r303-one-dimensional-mc クレートも common に依存しないスタンドアローン実装です。XorShift64 は 3.2 と同じものをローカルに持ちます。

PDF と逆関数法によるサンプリング関数は次の通りです。

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()
}

推定関数はそれぞれ $f(x)/p(x)$ を $N$ 回加算して平均を取ります。

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
}

1 回の実行では乱数の偶然に左右されるため，summarize_trials で同じ条件を複数回繰り返して平均絶対誤差を計算します。uniform と importance それぞれに独立した RNG インスタンスを使うことで，一方の試行結果が他方に影響しません。

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 は 2 段階の出力を構成します。まず $N={10}^6$ の単発推定で 2 手法の精度を比較し，続いて $N=100\sim 10000$ の小サンプルで複数試行の平均絶対誤差を CSV 形式で記録します。

C++ と Rust の違い

C++ 原典では rand() によるグローバルな乱数状態を使います。

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

Rust では rng を引数として明示的に渡します。summarize_trials が rng_u と rng_i を独立したインスタンスで保持することで，2 手法の試行系列が互いに干渉せず，再現性が保証されます。

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
}

ブラウザ上での実行結果

method,estimate,abs_error の行では $N={10}^6$ の単発推定で 2 手法を比較します。samples,trials,uniform,importance_x_over_2 の行では，左 2 列が試行条件，右 2 列が平均絶対誤差です。小サンプルほど差が大きく現れ，通常は importance_x_over_2 が小さくなります。

まとめ

• ${\hat{I}}_N=\frac{1}{N}\sum f(x_i)/p(x_i)$ というモンテカルロ推定量の一般形を導いた。$p$ が何であれ不偏推定量であり，$p$ の選び方が分散に影響する。
• 一様サンプリング $p(x)=1/2$ で ${\int }_0^2{x}^{2}dx$ を推定した。
• $p(x)=x/2$ による重点サンプリングは逆関数法 $x=\sqrt{4u}$ でサンプルでき，同じサンプル数で分散が下がる。
• summarize_trials による複数試行の平均絶対誤差で，1 回の試行では見えにくい分散の差を実測した。

次章では，この $f(x)/p(x)$ という構造を球面積分へ拡張します。