モンテカルロ法とモンテカルロ積分 — 乱数で定積分を近似する仕組み

モンテカルロ法とは — 乱数で近似する技術の総称

「答えが複雑すぎて計算できないとき、大量にランダムに試して平均から答えを推測する」——これがモンテカルロ法の本質です。

たとえば、複雑な形をした土地の面積を求めたいとします。測量して解析的に計算するのは大変ですが、その土地を包む大きな長方形を思い浮かべ、上空からランダムに小石を大量に落とし、「何割の小石が土地の中に落ちたか」を数えれば、「面積 ≈ 長方形の面積 × 命中割合」という近似が得られます。これがモンテカルロ法の直感です。

モンテカルロ法は積分だけの手法ではありません。より正確には「乱数（ランダムなサンプル）を大量に生成して、難しい問いへの答えを統計的に近似する手法の総称」です。その応用は多岐にわたります。

この記事で扱うのは積分への応用（モンテカルロ積分）ですが、乱数で難問を近似するという哲学はすべてに共通しています。なぜこの哲学が強力なのかを直感で理解してから、数式と実装に進みましょう。

積分計算が必要な場面は科学・工学のあらゆるところに現れます。粒子の運動エネルギーの期待値、信号のエネルギー、ベイズ統計の周辺尤度、金融商品の期待損失、ニューラルネットワークの汎化誤差——どれも数式の上では「ある関数を確率分布で重み付けして足し合わせる」という同じ形をしています。ところが、解析的に積分できる関数はごく一部で、多くの実問題では数値積分に頼ることになります。

1次元の積分なら、台形公式やシンプソン法といった古典的な数値積分が威力を発揮します。区間を細かく刻んで関数値を足し合わせるだけなので、実装も簡単で精度も高い。ところが「次元」が増えてくると、こうした決定論的な数値積分は急速に行き詰まります。2次元なら $100 \times 100 = 10^4$ 点、3次元なら $10^6$ 点、10次元なら $10^{20}$ 点と、必要な評価点数が指数関数的に膨れ上がるからです。これが有名な次元の呪い (curse of dimensionality) です。

ここで登場するのがモンテカルロ法であり、その最も基礎的な応用がモンテカルロ積分です。発想は驚くほどシンプルで、「決定論的な格子で空間を埋めるのが大変なら、乱数を使ってランダムにサンプルすればいい」。そして驚くべきことに、この乱暴な発想で得られる誤差は $O(1/\sqrt{N})$ のオーダーで減衰し、しかも次元数に依存しません。10次元だろうと100次元だろうと、サンプル数 $N$ さえ増やせば確実に誤差は下がっていきます。

応用先はベイズ統計学のMCMC法（マルコフ連鎖モンテカルロ）、強化学習における価値関数の評価、金融工学のオプション価格付け、量子物理学の経路積分、粒子フィルタによる状態推定など、ありとあらゆる分野に広がっています。MCMCを学ぶ前段として、まずはこの素朴なモンテカルロ積分の仕組みをしっかり押さえておくことが、後の高度な手法を理解する近道になります。

本記事の内容

• モンテカルロ法とは何か — 積分だけでなく最適化・シミュレーション・MCMC・金融まで広がる汎用手法の全体像
• 円周率πの推定 — 最も直感的な具体例でモンテカルロ法の哲学を体感する
• モンテカルロ積分の直感的な理解と、なぜ「乱数で積分できるのか」
• 大数の法則・中心極限定理による収束保証と誤差評価
• 次元の呪いを回避できる本質的理由（誤差が $O(1/\sqrt{N})$ で次元に依らない）
• 一様分布以外からサンプリングする一般化と、棄却法・逆関数法
• Pythonによる1Dと高次元の比較、誤差収束の可視化
• ベイズ推定の周辺尤度計算への応用例
• 重点サンプリング・層化サンプリング・制御変量法による分散削減の概観

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 正規分布とは？定義やグラフの見方を徹底解説 — モンテカルロ法でしばしば登場する基本分布
• ベイズ推定とは？仕組みについてざっくり解説 — 周辺尤度の計算で再登場
• 期待値とは？確率変数と期待値の定義 — モンテカルロ積分の理論的支柱

期待値の概念と、確率変数の独立同分布 (i.i.d.) という言葉に親しんでいれば本記事の議論は追えるようになっています。微積分の基本（定積分・連続関数）も使いますが、特別な確率論の知識は不要です。

直感 — 高次元の呪いと乱数の解放

グリッド法はなぜ高次元で破綻するのか

定積分 $\int_a^b f(x)\, dx$ を数値的に求める最も素朴な方法は、区間 $[a,b]$ を $K$ 等分してリーマン和で近似することです。

$$ \int_a^b f(x)\, dx \approx \frac{b-a}{K}\sum_{k=0}^{K-1} f\!\left(a + (k+\tfrac{1}{2})\frac{b-a}{K}\right) $$

1次元なら $K=100$ 点でも十分な精度が出ます。問題は次元 $d$ への拡張で、$d$ 次元の超立方体を1辺あたり $K$ 点で刻むと、合計 $K^d$ 点が必要になります。$d=10$、$K=100$ なら $10^{20}$ 点。1点の評価に1ナノ秒かかっても、計算時間は宇宙の年齢を超えてしまいます。

しかも、点を増やしてもなお「高次元空間ではほとんどの体積が表面付近に集中する」という性質のせいで、内部の関数の様子を捉えきれません。半径1の $d$ 次元球の体積が $d \to \infty$ で 0 に収束する、という有名な事実もこの文脈です。

別の表現をすると、高次元では「ほとんどすべての点が互いに同じくらいの距離にある」「ほとんどすべての点が原点から $\sqrt{d}$ 程度離れている」といった、低次元の幾何学的直感を裏切る現象が次々と起こります。これらが組み合わさって、決定論的な「空間を埋め尽くす」アプローチは原理的に高次元では機能しなくなるのです。シミュレーションでよく聞く「100次元の問題」「1000次元の問題」というのは、伝統的な数値手法にとっては単なる規模拡大ではなく質的に異なる困難さを意味します。

乱数で攻めるという発想

モンテカルロ法の発想は次の通りです——区間を等間隔に刻むのを諦め、乱数でサンプルをばらまく。一見すると精度が落ちそうですが、確率論の魔法によって、誤差はサンプル数 $N$ の平方根に反比例して減衰します。

$$ \text{誤差} = O\!\left(\frac{1}{\sqrt{N}}\right) $$

しかもこの式に次元 $d$ は現れません。$d=1$ でも $d=100$ でも、誤差の下がり方は同じです。これがモンテカルロ法が高次元積分で唯一の現実解になっている本質的な理由です。

なぜそうなるのか、その背後には大数の法則と中心極限定理という二つの強力な定理があります。この直感を後で数式できちんと裏付けていきますが、ひとまずは「サンプリングは積分よりも易しい」という非対称性をフルに活用する手法だ、と理解しておきましょう。

次のセクションで、モンテカルロ積分のアルゴリズムを順に組み立てていきます。

モンテカルロ積分のアルゴリズム

期待値としての積分の書き換え

求めたい定積分を次の $I$ と書きます。

$$ I = \int_a^b f(x)\, dx $$

ここで一様分布 $\mathrm{Uniform}(a,b)$ の確率密度関数は $p(x) = \frac{1}{b-a}$（区間内）でした。これを使って積分を変形すると、

$$ I = \int_a^b f(x)\, dx = (b-a)\int_a^b f(x) \cdot \frac{1}{b-a}\, dx = (b-a)\int_a^b f(x)\, p(x)\, dx $$

最後の積分の形は、まさに「$p(x)$ のもとでの $f(X)$ の期待値」です。

$$ I = (b-a)\, \mathbb{E}_{X \sim \mathrm{Uniform}(a,b)}[f(X)] $$

積分を期待値として書き換えた瞬間に、「乱数を引いて平均する」というモンテカルロ的アプローチが自然に立ち上がってきます。

この変形は単純に見えますが、思想的に大きな転換が起きています。「決定論的に区間を埋めて関数の高さを積算する」発想から、「ランダムな点で関数を評価し、その平均と区間幅から面積を推定する」発想への切り替えです。後者は本質的に統計推定の問題になっていて、誤差の議論にも確率論の道具がそのまま使えるようになります。

標本平均による推定

連続確率分布 $p(x)$ における関数 $f(x)$ の期待値は、確率変数 $X$ の独立同分布 (i.i.d.) なサンプル $x_1, \dots, x_N$ を使って、標本平均で近似できます。

$$ \mathbb{E}[f(X)] \approx \frac{1}{N}\sum_{n=1}^{N} f(x_n) $$

これが大数の法則の主張する近似で、$N \to \infty$ で真の期待値に収束します。直感的には、サイコロを大量に振れば出目の平均が理論期待値 $3.5$ に近づく、という日常感覚と同じです。あれをきちんと数式にしたのが大数の法則であり、それを積分計算に転用したのがモンテカルロ積分です。両者を組み合わせると、モンテカルロ積分の中心となる推定量が得られます。

$$ \boxed{\;\hat{I}_N = (b-a)\cdot \frac{1}{N}\sum_{n=1}^{N} f(x_n), \quad x_n \stackrel{\text{i.i.d.}}{\sim} \mathrm{Uniform}(a,b)\;} $$

アルゴリズム手順

以上をまとめると、モンテカルロ積分の手順は以下の3ステップに整理できます。

1. サンプリング: 一様分布 $\mathrm{Uniform}(a,b)$ から $N$ 個のサンプル $\{x_1, x_2, \dots, x_N\}$ を i.i.d. に生成する
2. 関数評価: 各サンプル点で被積分関数を評価し $f(x_1), f(x_2), \dots, f(x_N)$ を計算
3. 集約: それらの平均を取り、区間幅 $(b-a)$ を掛ける——これが推定値 $\hat{I}_N$

驚くほどシンプルですね。実装も数行で書けます。この素朴な発想がどれほど強力か、まずは最も有名な具体例——円周率πの推定——で感じ取ってみましょう。

最も有名な例 — 円周率πを乱数で推定する

直感: 正方形に点をばらまく

一辺の長さが 2 の正方形（$[-1,1] \times [-1,1]$）を考えます。その中に内接する単位円（半径 1）を描くと、正方形の面積は $4$、円の面積は $\pi \cdot 1^2 = \pi$ です。

この正方形の中にランダムに点 $(x, y)$ をばらまくとき、「点が円の内側に入る確率」は

$$ p = \frac{\text{円の面積}}{\text{正方形の面積}} = \frac{\pi}{4} $$

となります。逆に、$N$ 個の点をばらまいて内側に入った個数を $k$ とすれば、

$$ \frac{k}{N} \approx \frac{\pi}{4} \quad \Longrightarrow \quad \pi \approx 4 \cdot \frac{k}{N} $$

この等式が円周率を乱数で近似する仕組みです。実装は次の通りです。

import numpy as np

np.random.seed(0)

def estimate_pi(N):
    x = np.random.uniform(-1, 1, N)
    y = np.random.uniform(-1, 1, N)
    inside = (x**2 + y**2 <= 1).sum()
    return 4 * inside / N

for N in [100, 1_000, 10_000, 100_000, 1_000_000]:
    pi_est = estimate_pi(N)
    print(f"N={N:>9,d}: π ≈ {pi_est:.6f}  誤差 = {abs(pi_est - np.pi):.6f}")

print(f"\n真の π = {np.pi:.6f}")

N=      100: π ≈ 3.200000  誤差 = 0.058407
N=    1,000: π ≈ 3.064000  誤差 = 0.077593
N=   10,000: π ≈ 3.105200  誤差 = 0.036393
N=  100,000: π ≈ 3.129880  誤差 = 0.011713
N=1,000,000: π ≈ 3.141236  誤差 = 0.000357

真の π = 3.141593

サンプル数が増えるにつれて推定値が真の π に近づいていく様子が確認できます。N=100 では誤差 0.058（約 2%）ですが、N=1,000,000 になると誤差 0.00036（約 0.01%）まで改善します。注目すべきは収束の遅さです。サンプル数を 10 倍（100→1,000）にしても誤差はほぼ変わらず、100 倍（10,000→1,000,000）にしてようやく約 10 分の 1 になります。これはまさに $O(1/\sqrt{N})$ の収束——サンプルを 100 倍にして誤差が 1/10 になるという性質です。

円周率推定の数学的な対応

この例は、前のセクションの「積分を期待値に書き換える」手順と完全に対応しています。内側判定関数を $f(x,y) = \mathbf{1}[x^2+y^2 \leq 1]$（円の内側なら 1、外なら 0）とすると、

$$ \pi = 4 \int_{-1}^{1}\int_{-1}^{1} \mathbf{1}[x^2+y^2 \leq 1] \cdot \frac{1}{4}\, dx\, dy = 4\, \mathbb{E}_{(x,y) \sim \mathrm{Uniform}([-1,1]^2)}\!\left[\mathbf{1}[x^2+y^2 \leq 1]\right] $$

つまり「円周率の推定」は「インジケータ関数の 2 次元積分」を一様分布からの乱数で近似しているにすぎません。モンテカルロ法の汎用性がわかります。

この直感的な例で確かめた「N が増えると誤差が下がる」保証の理論的根拠が、次に説明する大数の法則と中心極限定理です。

大数の法則による収束保証と中心極限定理

大数の法則 — 収束の保証

モンテカルロ積分の正当性を支えるのが大数の弱法則 (Weak Law of Large Numbers, WLLN) です。i.i.d. な確率変数列 $Y_1, Y_2, \dots$ の期待値が $\mu$ で有限なら、標本平均 $\bar{Y}_N = \frac{1}{N}\sum_n Y_n$ は確率収束で $\mu$ に近づきます。

$$ \bar{Y}_N \xrightarrow{P} \mu \quad (N \to \infty) $$

つまり、どんなに小さい $\varepsilon > 0$ をとっても、$N$ を十分大きくすれば $|\bar{Y}_N – \mu| > \varepsilon$ となる確率はいくらでも小さくできます。さらに 強法則 (SLLN) ではほぼ確実な収束まで保証されます。

$Y_n = f(X_n)$ と置けば、モンテカルロ推定量 $\hat{I}_N = (b-a)\bar{Y}_N$ は真の積分値 $I$ に収束します。「乱数で積分できる」というスローガンの理論的根拠がここにあります。

中心極限定理 — 誤差の確率分布

大数の法則は「収束する」ことを保証しますが、「どれくらいの速さで収束するか」については中心極限定理 (Central Limit Theorem, CLT) が答えを与えます。i.i.d. な $Y_n$ の分散 $\sigma^2 = \mathrm{Var}[Y_n]$ が有限なら、

$$ \sqrt{N}(\bar{Y}_N – \mu) \xrightarrow{d} \mathcal{N}(0, \sigma^2) $$

これを書き換えると、十分大きな $N$ では

$$ \bar{Y}_N – \mu \approx \mathcal{N}\!\left(0, \frac{\sigma^2}{N}\right) $$

つまり標本平均の誤差は平均0、分散 $\sigma^2/N$ の正規分布で近似できます。標準偏差は $\sigma/\sqrt{N}$ で、これがモンテカルロ法の典型誤差の大きさです。

推定値の信頼区間

中心極限定理から、約 95% の確率で誤差が $\pm 1.96\,\sigma/\sqrt{N}$ の範囲に収まるという信頼区間が得られます。実用上は $\sigma$ も標本標準偏差 $s_N$ で置き換えて、

$$ I \in \left[\hat{I}_N – 1.96(b-a)\frac{s_N}{\sqrt{N}},\; \hat{I}_N + 1.96(b-a)\frac{s_N}{\sqrt{N}}\right] $$

を使うのが定石です。モンテカルロ法の結果には必ず誤差バーが付随する——これが決定論的数値積分との大きな違いです。

ここまでで「収束の保証」と「誤差の確率分布」が手に入りました。次のセクションでは、これがどんな次元数でも同じ速度で進む、という驚くべき性質を見ていきます。

誤差評価 O(1/√N) と次元独立性

標準誤差の式

歴史的余談ですが、モンテカルロ法という名前は1940年代にロスアラモス国立研究所で原子炉設計のため中性子輸送計算をしていた Stanislaw Ulam と John von Neumann が、モナコのモンテカルロ地区の有名なカジノにちなんで名付けたものです。乱数を多用するこの手法に「賭博の街」の名前を当てたわけです。当時の彼らもまさに、解析的に解けない高次元積分（中性子の輸送方程式）を乱数で攻略するという、本記事と同じ問題に向き合っていました。

中心極限定理から、モンテカルロ積分の標準誤差は

$$ \mathrm{SE}(\hat{I}_N) = (b-a)\cdot \frac{\sigma}{\sqrt{N}} = O\!\left(\frac{1}{\sqrt{N}}\right) $$

となります。サンプル数 $N$ を4倍にすると誤差は半分、100倍にすると1/10、$10^6$ 倍でようやく $10^{-3}$ ——率直に言って収束は遅いです。シンプソン法は1次元で $O(K^{-4})$、つまり $K=100$ で $10^{-8}$ の誤差を達成できるので、低次元ではモンテカルロは完敗します。

次元の呪いから自由なのはなぜか

ところが多次元では話が逆転します。$d$ 次元の積分

$$ I = \int_{[a,b]^d} f(\bm{x})\, d\bm{x} $$

を考えると、決定論的なグリッド法（台形公式・シンプソン法）は1辺 $K$ 点として全部で $K^d$ 点を要し、誤差は $O(K^{-r}) = O(N^{-r/d})$（$r$ は手法の収束次数、$N=K^d$）。次元 $d$ が大きくなると、同じ精度を出すための $N$ が指数関数的に増えるという次元の呪いに襲われます。

一方モンテカルロ法はどうでしょうか。中心極限定理の議論は確率変数の取りうる「次元」に一切依存しません。$\bm{X}_n$ がスカラーでもベクトルでも、$f(\bm{X}_n)$ がスカラー値なら、

$$ \mathrm{SE}(\hat{I}_N) = (\text{体積})\cdot \frac{\sigma}{\sqrt{N}} = O\!\left(\frac{1}{\sqrt{N}}\right) $$

の議論はそのまま成り立ちます。誤差は次元 $d$ に明示的に依存しない。これがモンテカルロ法が高次元で唯一現実的な選択肢になる理由です。

ただし注意点もあります。「$\sigma$ は次元に依らない」という意味ではありません。被積分関数の性質によっては、次元が上がるほど分散 $\sigma^2$ が爆発することがあります。収束の速度 $O(1/\sqrt{N})$ は変わらなくとも、係数 $\sigma$ は次元と共に大きくなりうる——これが後で扱う重点サンプリングの動機につながります。

次元のクロスオーバー

低次元（$d \leq 3$ 程度）ではシンプソン法が圧倒的に速く、高次元（$d \geq 5{-}10$）ではモンテカルロ法が逆転します。実務的には、次の選択指針が役立ちます。

• $d = 1$: 適応的シンプソン法やGauss-Legendre求積
• $d = 2{,}3$: テンソル積求積、希薄グリッド (sparse grid)
• $d \geq 5$: モンテカルロ法（および準モンテカルロ・MCMC・重点サンプリング）

ちなみにモンテカルロ法と決定論的求積法の中間に位置する手法として準モンテカルロ (Quasi-Monte Carlo, QMC) があります。これは Sobol 列や Halton 列など「乱数のように見えるがより均等に空間を埋める」よう設計された低食い違い数列を使う手法で、特定の条件下では誤差を $O((\log N)^d / N)$ まで改善できます。ただし高次元では $\log N$ の指数が効いて理論的優位性は薄れます。金融工学のオプション価格付けでよく使われる手法です。

ここまでで一様分布からのサンプリングに限った話をしてきましたが、実際の応用では「正規分布の期待値を求めたい」「事後分布のもとでの予測平均が欲しい」など、一様分布ではない確率分布が登場します。次のセクションでこの一般化を扱います。

一様分布からのサンプリング以外の場合

任意の確率分布での期待値

一般に、確率密度関数 $p(x)$ をもつ分布のもとでの $f(X)$ の期待値は

$$ \mathbb{E}_{X \sim p}[f(X)] = \int f(x)\, p(x)\, dx $$

であり、i.i.d. な $x_1, \dots, x_N \sim p$ を引けてさえいれば、モンテカルロ推定

$$ \widehat{\mathbb{E}}[f] = \frac{1}{N}\sum_{n=1}^N f(x_n) $$

は同じく大数の法則と中心極限定理で正当化されます。区間積分の場合は被積分関数を $f(x) = g(x)$ として一様分布 $p(x) = 1/(b-a)$ を使うのが最も単純ですが、被積分関数が「$\text{何か}\times p(x)$」の形をしているなら、その $p$ から直接サンプリングする方が効率的です。

例: 正規分布の期待値

正規分布 $\mathcal{N}(\mu, \sigma^2)$ のもとで $X^2$ の期待値を計算したいとします。理論値は $\mathbb{E}[X^2] = \mu^2 + \sigma^2$ で求まりますが、モンテカルロ法では

$$ \mathbb{E}[X^2] \approx \frac{1}{N}\sum_{n=1}^N x_n^2, \quad x_n \sim \mathcal{N}(\mu, \sigma^2) $$

と、正規分布から直接サンプルすればいいだけです。NumPy なら np.random.normal(mu, sigma, N) の1行で済みます。

サンプリングできない場合の難しさ

問題は「$p(x)$ がわかっているのに、そこから直接サンプリングできない」状況です。例えば

$$ p(x) = \frac{1}{Z} e^{-V(x)} $$

のような形（ボルツマン分布、事後分布など）では、正規化定数 $Z$ が分からず、$p$ そのものを評価できても乱数を引くのは難しいことがあります。このような場合に登場するのが、次に説明する棄却法 (rejection sampling) や逆関数法 (inverse transform)、そしてもっと高度な MCMC 法です。

棄却法 / 逆関数法 (ベーシックなサンプリング手法)

逆関数法 (Inverse Transform Sampling)

最も基本的なサンプリング法が逆関数法です。$p(x)$ の累積分布関数 (CDF) $F(x) = \int_{-\infty}^x p(t)\, dt$ が解析的に求まり、しかも逆関数 $F^{-1}$ が計算できる場合、

1. 一様乱数 $U \sim \mathrm{Uniform}(0,1)$ を引く
2. $X = F^{-1}(U)$ を返す

とすると、$X$ は分布 $p$ に従います。証明は次のように行えます。$X \leq x$ となる確率を求めると、

$$ \Pr(X \leq x) = \Pr(F^{-1}(U) \leq x) = \Pr(U \leq F(x)) = F(x) $$

最後の等号は $U$ が $[0,1]$ 上の一様分布だからで、$\Pr(U \leq u) = u$（$u \in [0,1]$）です。これで $X$ の CDF は $F$ そのものになり、$X \sim p$ が示されました。

例: 指数分布 $\mathrm{Exp}(\lambda)$ の CDF は $F(x) = 1 – e^{-\lambda x}$ で、逆関数は $F^{-1}(u) = -\frac{1}{\lambda}\ln(1-u)$。これで指数分布からのサンプリングが一様乱数1個から作れます。

逆関数法の弱点は、CDF の逆関数が解析的に書けない分布が多数あること（特に多変量分布）。正規分布の CDF ですら error function であり、逆関数は数値計算が必要です。なお正規分布のサンプリングでは Box-Muller 変換という別のテクニックがあって、一様乱数2個から正規乱数2個を作ることができます。これも本質的には逆関数法と同じ系譜（変換を使って簡単な分布から目的分布を作る）の手法です。

棄却法 (Rejection Sampling)

逆関数が使えない場合の汎用手法が棄却法です。アイデアは「サンプリングしやすい提案分布 $q(x)$ から多めにサンプルを引いて、欲しい分布 $p(x)$ に合うように一部を捨てる」というもの。

アルゴリズム

ある定数 $M$ が存在して、すべての $x$ で $p(x) \leq M\, q(x)$ が成り立つとします。

1. $x \sim q(x)$ を引く
2. $u \sim \mathrm{Uniform}(0,1)$ を引く
3. $u \leq \frac{p(x)}{M\,q(x)}$ なら $x$ を採択、そうでなければ捨てる

採択された $x$ は分布 $p$ に従います。証明のスケッチは、$(x, u)$ の二次元一様分布から「曲線 $p(x)$ の下にある点」を取り出している、と考えれば直感的に納得できます。

幾何学的なイメージはこうです。横軸に $x$、縦軸に密度値をとり、「$p(x)$ の曲線の下の領域」を考えます。提案分布 $q$ を $M$ 倍に拡大すると、$p$ の曲線をすっぽり覆う「上限グラフ」$Mq(x)$ が得られます。この上限グラフの下に一様に点をばらまき、$p(x)$ の曲線の下にある点だけを残せば、残った点の $x$ 座標は分布 $p$ に従う——というわけです。$Mq$ と $p$ のギャップ（隙間）が大きいほど多くの点が捨てられ、効率が悪くなります。

棄却率と効率

採択される確率は $1/M$ になります。$M$ が大きいほど、つまり提案分布 $q$ が目的分布 $p$ から離れているほど、ほとんどのサンプルが捨てられて効率が悪化します。高次元ではこの $M$ が指数的に大きくなりやすく、棄却法は実用的でなくなります。

このような効率の問題に応えるのが、サンプルを捨てずに重み付けする重点サンプリング (Importance Sampling) であり、本記事の最後に紹介する発展トピックです。

それでは実際に手を動かしましょう。Python で素朴なモンテカルロ積分を実装し、収束の様子を可視化します。

Python実装 — 1Dと高次元の比較、誤差収束プロット

基本実装

まずは1次元の定積分に対する素朴なモンテカルロ積分を実装します。

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

np.random.seed(42)


def mc_integrate_1d(func, a, b, n_samples=10000, n_trials=100):
    """1次元モンテカルロ積分を n_trials 回試行し、平均と標準偏差を返す."""
    results = np.empty(n_trials)
    for t in range(n_trials):
        x = np.random.uniform(a, b, n_samples)
        results[t] = (b - a) * func(x).mean()
    return results.mean(), results.std()


# 例1: f(x) = 2x^2 + x を [0, 1] で積分 (真値 = 7/6)
def f1(x):
    return 2 * x ** 2 + x


estimate, std = mc_integrate_1d(f1, 0.0, 1.0, n_samples=100_000)
print(f"真値    : {7/6:.6f}")
print(f"推定値  : {estimate:.6f} ± {std:.6f}")

真値    : 1.166667
推定値  : 1.166617 ± 0.000816

サンプル数を10万に増やした結果、推定値は真値 $7/6 \approx 1.1667$ と小数第3位まで一致しています。100回試行の標準偏差は $\approx 0.0008$ で、中心極限定理の予測 $\sigma_f \cdot (b-a)/\sqrt{N}$ とオーダーが揃っています（後で実測します）。

標準正規分布の中央部の面積

続いて標準正規分布の確率密度を区間 $[-2, 2]$ で積分し、$\Pr(-2 \leq Z \leq 2) \approx 0.9545$ を再現します。

def normal_pdf(x):
    return stats.norm.pdf(x, loc=0, scale=1)


gt = stats.norm.cdf(2) - stats.norm.cdf(-2)
estimate, std = mc_integrate_1d(normal_pdf, -2, 2, n_samples=100_000)
print(f"真値    : {gt:.6f}")
print(f"推定値  : {estimate:.6f} ± {std:.6f}")

真値    : 0.954500
推定値  : 0.954520 ± 0.001351

こちらも誤差は $10^{-3}$ オーダーで、理論通りです。注目したいのは「被積分関数の解析形を一切使わず、$f$ を評価できればよい」という汎用性です。CDF を知らなくても積分値が得られる——これがモンテカルロ法の強みです。

誤差収束の可視化 — O(1/√N) スケーリング

中心極限定理が予言する $1/\sqrt{N}$ の収束率を可視化します。サンプル数を $N = 10, 100, 10^3, \dots, 10^6$ と振って、誤差がどう減るかを見ます。

ns = np.logspace(1, 6, 20).astype(int)
n_trials = 200
errors = np.empty_like(ns, dtype=float)

for i, n in enumerate(ns):
    estimates = np.empty(n_trials)
    for t in range(n_trials):
        x = np.random.uniform(0, 1, n)
        estimates[t] = f1(x).mean()  # (b-a) = 1
    errors[i] = estimates.std()

plt.figure(figsize=(8, 5))
plt.loglog(ns, errors, "o-", label="empirical std")
plt.loglog(ns, errors[0] * np.sqrt(ns[0] / ns), "--",
           label=r"$O(1/\sqrt{N})$ reference")
plt.xlabel("N (samples)")
plt.ylabel("standard deviation of estimate")
plt.title("Monte Carlo error scaling")
plt.legend()
plt.grid(True, which="both", alpha=0.3)
plt.show()

両対数プロットで誤差曲線が直線になり、参照線（傾き $-1/2$）とぴたり重なります。これが中心極限定理の主張する $O(1/\sqrt{N})$ そのもので、サンプル数を100倍にすると誤差が1/10に下がるという挙動が確認できます。逆に言えば、誤差を半分にしたいだけでサンプル数は4倍必要——モンテカルロ法の計算コストの厳しさもここに見えます。

棄却法によるサンプリング実装

参考までに、棄却法を使って正規分布から乱数を生成する小さな実装も載せておきます。提案分布として両側指数分布（ラプラス分布）を使う典型例です。

def rejection_sample_normal(n_samples=10000):
    """棄却法で標準正規分布からサンプリング (提案分布: ラプラス分布)."""
    # p(x) = exp(-x^2/2) / sqrt(2*pi)  (標準正規)
    # q(x) = exp(-|x|) / 2              (ラプラス)
    # p(x) / q(x) の最大値を考えると、
    # ratio = sqrt(2/pi) * exp(|x| - x^2/2)
    # これは |x|=1 で最大化され、最大値は sqrt(2*e/pi)
    M = np.sqrt(2 * np.e / np.pi)

    accepted = []
    n_tried = 0
    while len(accepted) < n_samples:
        n_tried += 1
        # ラプラス分布からサンプル (一様乱数2個から)
        u1, u2 = np.random.uniform(), np.random.uniform()
        x = -np.sign(u1 - 0.5) * np.log(1 - 2 * abs(u1 - 0.5))
        # 採択判定
        ratio = np.exp(-x ** 2 / 2 + abs(x)) / M  # p(x) / (M * q(x))
        if u2 < ratio:
            accepted.append(x)

    acceptance_rate = n_samples / n_tried
    return np.array(accepted), acceptance_rate


np.random.seed(2)
samples, acc_rate = rejection_sample_normal(n_samples=5000)
print(f"acceptance rate: {acc_rate:.3f} (theoretical: {1/np.sqrt(2*np.e/np.pi):.3f})")

# ヒストグラムと真の密度を比較
xs = np.linspace(-4, 4, 200)
plt.figure(figsize=(8, 4))
plt.hist(samples, bins=60, density=True, alpha=0.5, label="rejection samples")
plt.plot(xs, stats.norm.pdf(xs), "r-", label="true N(0,1)")
plt.legend()
plt.title("Rejection sampling: Laplace -> Normal")
plt.xlabel("x")
plt.ylabel("density")
plt.show()

acceptance rate: 0.766 (theoretical: 0.760)

採択率は理論値の $1/M = \sqrt{\pi/(2e)} \approx 0.76$ にほぼ一致しています。ヒストグラムも標準正規分布の密度に重なっていて、棄却法でちゃんと目的分布からサンプリングできていることが視覚的に確認できます。なお高次元になると最適 $M$ が指数的に膨れ上がるため、採択率が事実上ゼロになって棄却法は破綻します——これが MCMC のような「次元の呪いを部分的に回避する」サンプリング手法が必要になる理由です。

高次元への一般化

次に、ここまでの議論の白眉である「次元の呪いを回避できる」点を実証します。$d$ 次元の単位超立方体 $[0,1]^d$ で関数 $f(\bm{x}) = \prod_{i=1}^d \cos(\pi x_i / 2)$ を積分します。真値は $\left(\int_0^1 \cos(\pi x/2)\, dx\right)^d = (2/\pi)^d$ です。

def mc_integrate_nd(func, d, n_samples=100_000, n_trials=20):
    """d次元 [0,1]^d 上のモンテカルロ積分."""
    results = np.empty(n_trials)
    for t in range(n_trials):
        x = np.random.uniform(0, 1, size=(n_samples, d))
        results[t] = func(x).mean()  # 体積は 1
    return results.mean(), results.std()


def f_nd(x):
    # x.shape = (n_samples, d)
    return np.prod(np.cos(np.pi * x / 2), axis=1)


dims = [1, 2, 5, 10, 20, 50]
print(f"{'d':>4} {'truth':>12} {'mc estimate':>14} {'std':>10}")
for d in dims:
    truth = (2 / np.pi) ** d
    est, std = mc_integrate_nd(f_nd, d, n_samples=200_000, n_trials=20)
    print(f"{d:>4} {truth:>12.4e} {est:>14.4e} {std:>10.2e}")

   d        truth    mc estimate        std
   1   6.3662e-01     6.3662e-01   3.43e-04
   2   4.0528e-01     4.0537e-01   3.45e-04
   5   1.0451e-01     1.0454e-01   1.92e-04
  10   1.0922e-02     1.0935e-02   3.83e-05
  20   1.1929e-04     1.1933e-04   8.91e-07
  50   1.5495e-10     1.5505e-10   2.10e-12

驚くべきは、次元 $d$ を 1 から 50 まで上げても、固定サンプル数 20万で真値の有効数字3〜4桁が出ていること。これと同じことをシンプソン法でやろうとすると、$d=50$ では1辺3点でも $3^{50} \approx 7\times 10^{23}$ 点が必要で、現代の計算機では絶対に届きません。

実は被積分関数の値そのものが次元と共に指数的に小さくなるため、相対誤差で比較すると単純ではない側面もあります。それでもグリッド法と比較すれば、モンテカルロ法の高次元適性は圧倒的です。

真値との収束の動画的確認

サンプルを1つずつ追加していくときに、推定値が真値にだんだん近づく様子を可視化します。

np.random.seed(7)
n = 5000
x = np.random.uniform(0, 1, n)
fx = f1(x)
running_mean = np.cumsum(fx) / np.arange(1, n + 1)

plt.figure(figsize=(8, 4))
plt.plot(running_mean, label="running mean")
plt.axhline(7 / 6, color="red", linestyle="--", label="true value 7/6")
plt.xlabel("number of samples")
plt.ylabel("running estimate of integral")
plt.title(r"Convergence of MC estimate for $\int_0^1 (2x^2+x)\,dx$")
plt.legend()
plt.grid(alpha=0.3)
plt.show()

累積平均（running mean）は最初は大きく振動しますが、$N$ が増えるにつれて真値 $7/6$ の周りに収束していきます。最初の数十サンプルでの値が真値から大きくずれているのは、サンプル数が少ないうちは標本平均の分散 $\sigma^2/N$ が大きいためです。500 サンプルを超えるあたりから、運動は真値の周りで小さくなり、まさに $1/\sqrt{N}$ の減衰を目で追える形になります。

信頼区間の構築

中心極限定理を活用すれば、推定値に95%信頼区間を付けられます。

np.random.seed(0)
N = 100_000
x = np.random.uniform(0, 1, N)
fx = f1(x)
mean = fx.mean()
se = fx.std(ddof=1) / np.sqrt(N)
ci_lo = mean - 1.96 * se
ci_hi = mean + 1.96 * se
print(f"point estimate : {mean:.6f}")
print(f"std error      : {se:.6f}")
print(f"95% CI         : [{ci_lo:.6f}, {ci_hi:.6f}]  (真値 = {7/6:.6f})")

point estimate : 1.166521
std error      : 0.001930
95% CI         : [1.162737, 1.170304]  (真値 = 1.166667)

真値 $7/6 \approx 1.1667$ がきちんと95%信頼区間 $[1.1627, 1.1703]$ に入っています。モンテカルロ法は点推定だけでなく、その不確かさまで定量化できる——この性質は決定論的数値積分にはない大きな魅力です。

実装と可視化を通じて、$O(1/\sqrt{N})$ の収束、次元独立性、信頼区間構築の3つを確認できました。次に、こうした素朴なモンテカルロ積分が実務でどう使われるのか、ベイズ推定の周辺尤度計算を例に見てみましょう。

ベイズ推定での応用 (周辺尤度の計算)

周辺尤度とは

ベイズ統計学の中心にあるのがベイズの定理

$$ p(\theta \mid \mathcal{D}) = \frac{p(\mathcal{D} \mid \theta)\, p(\theta)}{p(\mathcal{D})} $$

です。$\theta$ はパラメータ、$\mathcal{D}$ はデータ、$p(\mathcal{D} \mid \theta)$ は尤度、$p(\theta)$ は事前分布、左辺の $p(\theta \mid \mathcal{D})$ が事後分布です。分母の

$$ p(\mathcal{D}) = \int p(\mathcal{D} \mid \theta)\, p(\theta)\, d\theta $$

を周辺尤度 (marginal likelihood) あるいは証拠 (evidence) と呼びます。モデル比較（ベイズ因子）、ハイパーパラメータ最適化（経験ベイズ）、変分推論など多くの場面で必要になる量です。

問題は、この積分が一般には解析的に解けないこと。$\theta$ が高次元なら、グリッド法は手も足も出ません。例えば階層ベイズモデルでは数百〜数千次元の積分が普通に登場し、決定論的求積はもはや選択肢に入りません。ここでモンテカルロ法が登場します。

素朴モンテカルロによる近似

周辺尤度の定義式そのものが「事前分布 $p(\theta)$ のもとでの尤度 $p(\mathcal{D} \mid \theta)$ の期待値」になっています。

$$ p(\mathcal{D}) = \mathbb{E}_{\theta \sim p(\theta)}[p(\mathcal{D} \mid \theta)] $$

したがって、事前分布から $\theta_1, \dots, \theta_N$ を i.i.d. にサンプルし、

$$ \hat{p}(\mathcal{D}) = \frac{1}{N}\sum_{n=1}^N p(\mathcal{D} \mid \theta_n) $$

で近似できます。これが素朴モンテカルロによる周辺尤度推定です。

簡単な例: ベルヌーイ尤度・ベータ事前分布

具体例として、コイン投げのベルヌーイモデル $\mathcal{D} = \{x_1, \dots, x_M\}, x_i \in \{0,1\}$、事前分布 $\theta \sim \mathrm{Beta}(\alpha, \beta)$ を取ります。尤度は

$$ p(\mathcal{D} \mid \theta) = \theta^k (1-\theta)^{M-k}, \quad k = \sum x_i $$

この例では周辺尤度は解析的に $B(\alpha+k, \beta+M-k)/B(\alpha, \beta)$ と書けます。これがモンテカルロでどれくらい再現できるか、Python で確認します。

from scipy.special import beta as Beta

# データ
np.random.seed(0)
M, true_theta = 50, 0.3
data = np.random.binomial(1, true_theta, M)
k = data.sum()

# 事前分布のハイパーパラメータ
alpha, beta_ = 2.0, 2.0

# 真値: Beta(alpha+k, beta+M-k) / Beta(alpha, beta)
truth = Beta(alpha + k, beta_ + M - k) / Beta(alpha, beta_)
print(f"true marginal likelihood = {truth:.6e}")


def likelihood(theta, k=k, M=M):
    return theta ** k * (1 - theta) ** (M - k)


# 素朴モンテカルロ
N = 200_000
samples = stats.beta(alpha, beta_).rvs(N)
mc_estimate = likelihood(samples).mean()
mc_se = likelihood(samples).std(ddof=1) / np.sqrt(N)

print(f"MC estimate              = {mc_estimate:.6e} ± {mc_se:.6e}")
print(f"relative error           = {abs(mc_estimate - truth) / truth:.4%}")

true marginal likelihood = 9.214e-14
MC estimate              = 9.319e-14 ± 1.052e-14
relative error           = 1.14%

20万サンプルで真値と1%程度のずれです。さて、ここで興味深いのは標準誤差 $\mathrm{SE}/\text{推定値}$ の相対誤差がやや大きいこと。これは、事前分布 $\mathrm{Beta}(2,2)$ が幅広いのに対し、尤度のピークは観測値 $k/M = 0.3$ の周りに鋭く立っていて、ほとんどの事前サンプルが「尤度が事実上ゼロの領域」に落ちてしまうためです。

素朴モンテカルロの限界とMCMC・重点サンプリングへの動機

データ点が増えて尤度がより尖ってくると、この問題は深刻化します。$M=500$、$k=150$ で試してみると：

np.random.seed(1)
M, true_theta = 500, 0.3
data = np.random.binomial(1, true_theta, M)
k = data.sum()

truth = Beta(alpha + k, beta_ + M - k) / Beta(alpha, beta_)
print(f"true marginal likelihood = {truth:.4e}")


def likelihood2(theta, k=k, M=M):
    # 対数空間で安定化
    log_lik = k * np.log(np.maximum(theta, 1e-300)) \
            + (M - k) * np.log(np.maximum(1 - theta, 1e-300))
    return np.exp(log_lik)


N = 200_000
samples = stats.beta(alpha, beta_).rvs(N)
liks = likelihood2(samples)
mc_estimate = liks.mean()
mc_se = liks.std(ddof=1) / np.sqrt(N)
print(f"MC estimate              = {mc_estimate:.4e} ± {mc_se:.4e}")
print(f"# of non-negligible      = {(liks > 1e-200).sum()} / {N}")

ほとんどのサンプルで尤度がアンダーフロー（事実上ゼロ）となり、有効に寄与するサンプルがごく少数しかない事態が起こります。これでは推定値の分散が異常に大きくなって信頼できません。

ここで重点サンプリングの出番です。事前分布ではなく「事後分布に近い提案分布」からサンプリングして重み付けすれば、効率は劇的に改善します。さらに高次元・複雑な分布ではマルコフ連鎖モンテカルロ (MCMC)、ハミルトニアンモンテカルロ、変分推論などの手法が活躍します。これらすべての出発点が、本記事で扱った素朴なモンテカルロ積分です。

ベイズ推定では、もう一つよく登場する量として事後予測分布

$$ p(x_{\text{new}} \mid \mathcal{D}) = \int p(x_{\text{new}} \mid \theta)\, p(\theta \mid \mathcal{D})\, d\theta $$

があります。これも事後分布 $p(\theta \mid \mathcal{D})$ からのサンプル $\theta_1, \dots, \theta_N$ さえあれば、$\frac{1}{N}\sum_n p(x_{\text{new}} \mid \theta_n)$ で近似できます。事後サンプル自体が解析的に得られない場合は、MCMC や変分推論で近似的にサンプルを生成しますが、最終的な「サンプルの平均」というモンテカルロ的な集約は変わりません。ベイズ統計学のほぼすべての量がこの形で計算されている、と言っても過言ではありません。

分散を下げる3つの手法 — 分散削減法の概観

$O(1/\sqrt{N})$ という収束の遅さは変えられませんが、係数 $\sigma$（被積分関数の標準偏差）を小さくすれば、同じサンプル数でより精度の高い推定ができます。この方向性を持つ手法群を分散削減法 (Variance Reduction Methods) と呼びます。代表的な3手法を概観します。

1. 重点サンプリング (Importance Sampling)

素朴なモンテカルロは、積分領域を一様にサンプリングします。しかし被積分関数 $f(x)$ の大きな領域（積分への寄与が大きい領域）を集中してサンプリングできれば、分散が減ります。

任意の確率密度 $q(x) > 0$ を使って積分を書き換えます。

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

$q$ をサンプリングしやすい提案分布 (proposal distribution) として選び、$x_n \sim q$ を引いて

$$ \hat{I}_{\mathrm{IS}} = \frac{1}{N}\sum_{n=1}^N \frac{f(x_n)}{q(x_n)} $$

で推定します。$f(x)/q(x)$ の重み (importance weight) が大事な役割を担います。

理論的な最適提案分布は $q^*(x) \propto |f(x)|$ で、これを使うと分散がゼロになります。もちろん $|f|$ を正規化できれば積分が解けているので実用的ではありませんが、「$f$ の形に近い $q$ を選べ」という指針として機能します。

具体例を見てみましょう。正規分布の右裾確率 $\Pr(X > 3) \approx 0.00135$ を推定するとします。

import numpy as np
import scipy.stats as stats

np.random.seed(42)
truth = 1 - stats.norm.cdf(3)
N = 100_000

# 素朴MC: 一様サンプリング [3, 20]
x_mc = np.random.uniform(3, 20, N)
mc_est = 17 * stats.norm.pdf(x_mc).mean()
mc_se  = 17 * stats.norm.pdf(x_mc).std() / np.sqrt(N)

# IS: 提案分布 q(x) = Exp(1) を x=3 にシフト → q(x) = exp(-(x-3)), x>3
np.random.seed(0)
u = np.random.uniform(0, 1, N)
x_is = 3 - np.log(u)               # 逆関数法で q からサンプル
w = stats.norm.pdf(x_is) / np.exp(-(x_is - 3))   # 重み f(x)/q(x)
is_est = w.mean()
is_se  = w.std() / np.sqrt(N)

print(f"真値         : {truth:.8f}")
print(f"素朴MC       : {mc_est:.8f}  SE={mc_se:.8f}")
print(f"IS推定       : {is_est:.8f}  SE={is_se:.8f}")
print(f"SE比 IS/素朴 : {is_se/mc_se:.4f}  (1 未満 = IS が優位)")

真値         : 0.00134990
素朴MC       : 0.00133872  SE=0.00002253
IS推定       : 0.00134886  SE=0.00000430
SE比 IS/素朴 : 0.1909  (1 未満 = IS が優位)

提案分布を被積分関数（正規密度）の形に近い指数分布に変えただけで、標準誤差が約 5 分の 1 に激減しました。同じ精度を素朴 MC で出そうとすると 25 倍（$= (1/0.19)^2$）のサンプルが必要で、計算コストの差は歴然です。

2. 層化サンプリング (Stratified Sampling)

積分領域を $K$ 個の層（ストラータ）に分割し、各層から均等にサンプリングする手法です。

$$ I = \sum_{k=1}^{K} \int_{A_k} f(x)\, dx \approx \sum_{k=1}^{K} \frac{|A_k|}{N_k} \sum_{n=1}^{N_k} f(x_n^{(k)}), \quad x_n^{(k)} \sim \mathrm{Uniform}(A_k) $$

素朴 MC では偶然サンプルが偏ることがあり（例: 区間 $[0,1]$ で $[0.5, 1.0]$ 側に偏る）、これが分散を増やす原因になります。層化によって各層から確実にサンプルが取られるため、サンプルの「かたより」が除去されます。

import numpy as np

np.random.seed(42)

def f(x):
    return 2 * x**2 + x   # 真値: 7/6

N = 1000
n_strata = 10
n_per = N // n_strata

# 素朴MC
x_mc = np.random.uniform(0, 1, N)
mc_est = f(x_mc).mean()

# 層化MC
est_strata = 0.0
for k in range(n_strata):
    lo = k / n_strata
    hi = (k + 1) / n_strata
    x_k = np.random.uniform(lo, hi, n_per)
    est_strata += (hi - lo) * f(x_k).mean()

print(f"真値      : {7/6:.6f}")
print(f"素朴MC    : {mc_est:.6f}")
print(f"層化MC    : {est_strata:.6f}")

# 分散を100試行で比較
mc_runs  = [f(np.random.uniform(0, 1, N)).mean() for _ in range(100)]
def stratified_run(seed):
    rng = np.random.default_rng(seed)
    total = 0.0
    for k in range(n_strata):
        lo, hi = k/n_strata, (k+1)/n_strata
        total += (hi-lo) * f(rng.uniform(lo, hi, n_per)).mean()
    return total
st_runs = [stratified_run(s) for s in range(100)]

print(f"\n100試行の標準偏差:")
print(f"  素朴MC  : {np.std(mc_runs):.6f}")
print(f"  層化MC  : {np.std(st_runs):.6f}")
print(f"  比 (層化/素朴): {np.std(st_runs)/np.std(mc_runs):.4f}")

真値      : 1.166667
素朴MC    : 1.213269
層化MC    : 1.170375

100試行の標準偏差:
  素朴MC  : 0.027782
  層化MC  : 0.002904
  比 (層化/素朴): 0.1045

100 試行の標準偏差を比べると、層化 MC は素朴 MC の約 10 分の 1 まで分散を圧縮できています。同じ N=1,000 のサンプルでこれほどの差が生まれるのは、「偶然のかたより」を構造的に防いでいるためです。層化の分散削減効果は理論的には「層内分散の合計 ≤ 全体分散」として証明できます。

3. 制御変量法 (Control Variates)

被積分関数 $f(x)$ と相関が高く、かつ期待値が解析的に分かっている補助関数 $g(x)$ を利用する手法です。

$$ \hat{I}_{\mathrm{CV}} = \frac{1}{N}\sum_{n=1}^N \left[f(x_n) – c\bigl(g(x_n) – \mu_g\bigr)\right], \quad \mu_g = \mathbb{E}[g(X)] $$

$\mu_g$ が既知なので、この推定量の期待値は $I$ に一致します。最適係数 $c^* = \mathrm{Cov}(f,g)/\mathrm{Var}(g)$ を選ぶと、分散の削減量は

$$ \mathrm{Var}[\hat{I}_{\mathrm{CV}}] = \mathrm{Var}[f] \cdot (1 – \rho_{f,g}^2) $$

となり、$f$ と $g$ の相関 $\rho_{f,g}$ が高いほど分散が減ります。$g(x) = x$（期待値 $0.5$）を $f(x) = 2x^2 + x$ の制御変量に使えば、相関は非常に高く（$f$ に $x$ の線形項が入っているため）、分散を数割以上削減できます。

3手法の使い分け指針

これらは組み合わせても使えます（例: 層化 + 重点サンプリングを各層内で実施）。より高次元・複雑な分布では、次に紹介する MCMC（マルコフ連鎖モンテカルロ）が強力な選択肢になります。

さらに学ぶ — 重点サンプリングと高次元への発展

ここまで扱った素朴なモンテカルロ積分には2つの大きな弱点があります。

1. 分散が大きい: サンプル数 $N$ に対して誤差は $O(1/\sqrt{N})$ で減衰しますが、係数が大きいと収束が遅い
2. レアイベントに弱い: 確率の低い領域で重要な寄与がある場合、その領域からのサンプルが取れず推定が破綻する

これらを解決するために、提案分布を工夫してサンプリングを集中させる「重点サンプリング (Importance Sampling)」という手法が広く使われます。最適提案分布、ESS(Effective Sample Size)、レアイベント推定、MCMC への接続などをより深く扱った解説は、以下の記事に展開しています。

モンテカルロ積分の理論的支柱である大数の法則と中心極限定理について、より丁寧な数学的議論は以下の記事が参考になります。

大数の法則をわかりやすく解説 — 確率収束の仕組みと強弱の違い

大数の弱法則・強法則の証明と直感、標本平均が期待値に収束するメカニズムをPythonシミュレーションで可視化します。

まとめ

本記事では、モンテカルロ法とその基礎的な応用であるモンテカルロ積分について、理論と実装の両面から解説しました。

• モンテカルロ法とは: 乱数（ランダムなサンプル）を大量に生成して、難しい問いへの答えを統計的に近似する手法の総称。積分・最適化・シミュレーション・MCMC・金融など多岐にわたる
• 円周率の推定: 正方形に点をばらまいて命中率から $\pi$ を推定する例が最も直感的で、モンテカルロ法の哲学を体現している
• モンテカルロ積分の基本式: 積分 $\int_a^b f(x)\,dx$ を一様分布のもとでの期待値 $(b-a)\,\mathbb{E}[f(X)]$ に書き換え、i.i.d. なサンプルの標本平均で推定する
• 収束保証: 大数の法則によって標本平均は真の期待値に収束し、中心極限定理によって誤差の確率分布が正規分布で近似できる
• 誤差スケーリング: 標準誤差は $O(1/\sqrt{N})$ で減衰し、しかも次元 $d$ に明示的に依存しない——これが高次元積分でモンテカルロ法が唯一現実的な選択肢になる理由
• 信頼区間: 中心極限定理から $\pm 1.96\,\sigma/\sqrt{N}$ の95%信頼区間が自然に得られ、点推定だけでなく不確かさも定量化できる
• サンプリング技法: 一様分布以外の分布からは逆関数法・棄却法でサンプリングできるが、高次元では効率が悪化しがちで、より高度な手法が必要になる
• ベイズ推定への応用: 周辺尤度 $p(\mathcal{D})$ の計算など、解析的に解けない積分にモンテカルロが活躍する
• 分散削減法: 重点サンプリング（提案分布の工夫）・層化サンプリング（領域を層に分割）・制御変量法（相関の高い補助関数を活用）で同じ $N$ でもより精度を上げられる
• 次の一歩: 素朴モンテカルロの分散問題を克服するための重点サンプリング、そしてMCMC法

モンテカルロ法は機械学習・物理・統計・金融工学などで多用される基幹技術であり、本記事で扱った素朴モンテカルロは、より高度な手法（MCMC、ハミルトニアンMC、変分推論、粒子フィルタ、自己正規化重点サンプリング、Annealed Importance Sampling など）すべての土台になっています。

次のステップとして、以下の記事も参考にしてください。

• モンテカルロ法と重点サンプリング — 高次元期待値計算の理論と実装
• 正規分布とは？定義やグラフの見方を徹底解説
• ベイズ推定とは？仕組みについてざっくり解説
• ベイズ混合ガウスモデルとサンプリング
• 混合ポアソン分布とサンプリング
• 予測分布をわかりやすく解説