2. 簡単なモンテカルロプログラム

Ray Tracing: The Rest of Your Life (v3.2.3): 2 A Simple Monte Carlo Program / 3.2 簡単なモンテカルロプログラム

モンテカルロ法（MC 法）の本質的な特徴は「単純なプログラムで真の値を統計的に近似でき，実行時間が長いほど精度が上がる」という点です。本章ではその最小例として $\pi$ の推定を題材にし，第 3 編全体の理論基盤を作ります。

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

• 乱数サンプルを使ったモンテカルロ推定量の構成
• 収束速度（$O(1/\sqrt{N})$）の実測確認
• 層化サンプリングによる分散低減

π の推定

$[-1,1]\times [-1,1]$ の正方形内に一様な乱数点 $(x,y)$ を打ち，単位円 $x^2+y^2<1$ に入った割合を数えます。正方形の面積は $4$，単位円の面積は $\pi$ なので

$$\frac{\text{内側の点数}}{N}\approx \frac{\pi }{4}\Rightarrow \pi \approx 4\frac{\text{内側の点数}}{N}$$

となります。確率変数の期待値を標本平均で代替するという構造がそのまま現れており，これがモンテカルロ積分の最も単純な例です。

サンプル数 $N$ が大きくなるほど推定は真値に近づきますが，収束速度は $O(1/\sqrt{N})$ と遅いです。$N$ を $100$ 倍にしても精度は $10$ 倍にしかなりません。

Rust 実装

この章は ray tracing の幾何学とは独立した確率論の準備であるため，r302-monte-carlo クレートは common に依存しないスタンドアローン実装としています。

固定シードによる再現性を確保するため，Xorshift64 乱数生成器をローカルに実装します。common::utils にも同等の実装がありますが，本章を読む段階では common の内部実装には触れていないため，ここで全体を見せます。

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

推定関数は $N$ 回のサンプルを打ち，内側の点数から $\pi$ を推定します。

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 は 2 段階の出力を構成します。まず N=1000 の単発推定で 1 回の試行がどの程度の精度になるかを示し，続いて $N=10$ から ${10}^6$ までのチェックポイントで収束の様子を CSV 形式で記録します（第 3 のブロックは次節の層化サンプリング比較です）。

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
}

層化サンプリング

収束速度 $O(1/\sqrt{N})$ の遅さは，純粋な乱数がもつ「偶然の偏り」から来ています。層化サンプリング（stratified sampling）はサンプル空間を均等なセルに分割し，各セル内から 1 点だけサンプリングすることでこの偏りを排除します。

$[-1,1{]}^2$ を $\sqrt{N}\times \sqrt{N}$ のグリッドに分割し，$(i,j)$ 番目のセル内から一様乱数 $(u,v)\in [0,1{)}^2$ で 1 点を取ります。

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

すべてのセルに必ず 1 点が入るため，サンプルの被覆が均一になり分散が下がります。非層化では偶然に点が集中する場所や空白が生まれますが，層化ではそれが起きません。

generate_report の最後のブロックで，同じ総サンプル数（$\sqrt{N}=1000$，計 ${10}^6$ サンプル）を使った非層化と層化の比較を出力します。各 $(i,j)$ 反復で非層化サンプルと層化サンプルを 1 点ずつ取ることで，条件を揃えた比較になります。

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"
    ));

同じ ${10}^6$ サンプルでも，層化サンプリングのほうが誤差が小さくなる傾向があります。ただしこの優位性は次元が上がるほど薄れます（次元の呪い）。2 次元の $\sqrt{N}\times \sqrt{N}$ 分割は 3 次元なら $N^{1/3}\times N^{1/3}\times N^{1/3}$ になり，同じ総サンプル数で覆えるセル数が急減します。実用レイトレーサでは，層化はピクセルサンプリングのような低次元の部分にのみ適用し，多次元積分の内側には通常の乱数を使うのが一般的です。

C++ と Rust の違い

C++ 原典では rand() でグローバルな乱数を生成します。

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

Rust には rand() のようなグローバル乱数状態がありません。本章では common に依存しないため Xorshift64 をローカルに実装し，rng を引数として明示的に渡しています。関数の副作用が引数と戻り値だけに限定されるこのスタイルは，再現性の確保にも有利です。

出力フォーマットについては，C++ が std::cout << std::fixed << std::setprecision(12) を使うのに対して，Rust は 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
}

ブラウザ上での実行結果

収束チェックポイントでは abs_error が単調には減らず，乱数の偏りで前後しながらも全体として小さくなっていく様子が確認できます。層化サンプリング比較では，同じ ${10}^6$ サンプルを使っても層化のほうが誤差が小さくなる傾向があります（乱数次第で逆転することもあります）。

まとめ

• モンテカルロ推定量を $\pi$ 推定の最小例で実装した。収束速度は $O(1/\sqrt{N})$ であり，$N$ を $100$ 倍にして精度は $10$ 倍にしかならない。
• 収束チェックポイントにより，サンプル数と精度の関係を実測した。
• 層化サンプリングは同じサンプル数でより均一な被覆を実現し，2 次元の場合に分散を低減できる。次元が上がるほど恩恵は薄れる（次元の呪い）。

次章では，この $\pi$ 推定を出発点として，積分対象 $f(x)$ の期待値を確率変数の平均として捉える一般的なモンテカルロ積分の枠組みに進みます。