ある機器システムで稀な故障が今後30日以内に発生する確率を知りたい——とします。物理モデルは既にあり、部品の特性値や動作条件の不確かさはガウス分布で表せます。ところが、故障という事象はあまりにも稀で、典型的には $10^{-6}$ オーダー。素直に乱数を100万回振ってシミュレーションしても、故障は数回しか起きず、確率の推定値は荒っぽい上にときどき0になってしまいます。
似たような問題は他にもあります。金融機関が翌日の損失額が一定値を超える確率（Value at Risk, VaR）を見積もるとき、ベイズ統計学者がモデルの周辺尤度（証拠）を計算するとき、強化学習で「もしあのとき別のポリシーを取っていたら」を評価するとき——どれも、ある高次元空間における確率分布のもとで「ある量の期待値」を計算するという同じ構造の問題に行き当たります。
問題は、対象の積分が解析的に解けない上に、空間が高次元すぎてグリッドで刻むことすら絶望的だという点です。ここで活躍するのがモンテカルロ法であり、その分散を劇的に減らす重点サンプリング (Importance Sampling) です。重点サンプリングを上手に使えば、稀少故障確率の推定で1万倍の効率改善も珍しくありません。これは「100万サンプル必要だったところを100サンプルで済む」ことを意味します。
応用はそこにとどまりません。ベイズ推定における周辺尤度の計算、変分推論の証拠下界 (ELBO)、強化学習のオフポリシー評価、レアイベントの信頼性解析、物理学の経路積分量子化——いずれも背後で同じ重点サンプリングの考え方が動いています。
左図では素朴MCが目的分布 $p$ からサンプルを引くため、稀少事象の領域にほとんど届きません。右図では提案分布 $q$（稀少域に集中）からサンプルし、重み $w=p/q$ で補正することで、少ないサンプルで精度の高い推定が実現できます。重みの大きさが「$p$ に比べて $q$ がサンプルを引きすぎた分を割り引く」という補正の役割を担います。
本記事の内容
- モンテカルロ法の収束理論——分散 $\sigma^2/N$、中心極限定理、信頼区間、そして「次元の呪い」を回避するという驚くべき性質
- 重点サンプリング (IS) の定式化と分散低減の仕組み
- 最適提案分布 $q^*(x) \propto |f(x)|p(x)$ の導出と、有効サンプル数 (ESS) による品質評価
- レアイベント推定における劇的な効率改善
- Self-normalized IS、weight degeneracy、Annealed Importance Sampling などの実用上の罠と対策
- Pythonでの素朴MC vs IS の分散比較、$P(X>4)$ 推定の1000倍効率化、ESS の劣化可視化、MCMCへの接続
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- モンテカルロ法とモンテカルロ積分の基礎 — 本記事の前提となる素朴モンテカルロの仕組み
- 正規分布とは？定義やグラフの見方を徹底解説
- 指数分布をわかりやすく解説
- ベイズ推定とは？仕組みについてざっくり解説
- ベイズ混合ガウスモデルとサンプリング
- 予測分布をわかりやすく解説
直感 — 次元の呪いとモンテカルロの強さ
高次元の積分がなぜ難しいか、まず素朴な数値積分から考えてみましょう。1次元区間 $[0,1]$ を $K=100$ 等分してリーマン和で積分を近似すれば、$K=100$ 点で十分な精度が出ます。これを $d$ 次元の超立方体 $[0,1]^d$ に拡張すると、必要な点数は $K^d$ 個。$d=10$ なら $10^{20}$ 点、$d=100$ なら宇宙の原子数を超えます。これが次元の呪い (curse of dimensionality) です。
ところがモンテカルロ法の誤差は、後で示すように $O(1/\sqrt{N})$ で減ります——次元 $d$ に依らずに。$N=10^6$ 個のサンプルを振れば、1次元だろうが100次元だろうが、誤差はだいたい $10^{-3}$ のオーダーに落ち着きます。これがモンテカルロ法が高次元で「使い物になる」ほとんど唯一の方法である理由です。
別の見方をしてみます。100次元空間で「正規分布の典型的な値」を見ようとして、100次元のグリッドで埋めるのは不可能ですが、100次元の正規分布から1000サンプル引くのはほぼ瞬時です。サンプリングは積分よりも遥かに易しい——この非対称性をフルに使うのがモンテカルロ法の哲学です。
ただし、「サンプリングできる」と「効率よくサンプリングできる」は別物です。とくに、関心のある領域が分布の極端な裾にある場合（機器故障や異常接近のような稀少事象）、素朴にサンプリングしても肝心の領域にはほとんどサンプルが落ちません。ここから先で扱う重点サンプリングは、まさにこの「サンプルが届かない問題」を解決するテクニックです。
まずはモンテカルロ法そのものの理論的な収束を確認し、その後で重点サンプリングを導入していきます。
モンテカルロ法の理論と収束
期待値計算としての積分
モンテカルロ法の中心にあるのは、積分を期待値として見直すという発想の転換です。確率密度関数 $p(x)$ を持つ確率変数 $X$ について、関数 $f$ の期待値は
$$ I = \mathbb{E}_p[f(X)] = \int f(x)\, p(x)\, dx $$
と書けます。逆に、どんな積分 $\int g(x)\, dx$ も、$p(x)$ を一つ選んで $g(x) = f(x)p(x)$ となるように $f$ を取り直せば期待値として再解釈できます。たとえば $[0,1]$ 上の $\int_0^1 g(x)\, dx$ は、$p(x) = 1$（一様分布）として $\mathbb{E}_{\mathrm{Unif}[0,1]}[g(X)]$ と等価です。
積分を期待値に書き換える利点は、大数の法則が使えることです。$X_1, \dots, X_N$ を $p$ から独立に取った標本とすると、標本平均
$$ \hat{I}_N = \frac{1}{N}\sum_{i=1}^{N} f(X_i) $$
は $N \to \infty$ で真値 $I$ に確率収束します。これが素朴モンテカルロ推定量 (crude Monte Carlo) です。$N$ を増やすほど精度が上がる、というシンプルで力強いアイデアです。
不偏性と分散
$\hat{I}_N$ の良し悪しを正確に評価するため、期待値と分散を見てみます。期待値は線形性から
$$ \mathbb{E}[\hat{I}_N] = \frac{1}{N}\sum_{i=1}^{N} \mathbb{E}[f(X_i)] = \mathbb{E}_p[f(X)] = I $$
となり、$\hat{I}_N$ は不偏推定量です。分散は $X_i$ が独立であることを使うと
$$ \mathrm{Var}(\hat{I}_N) = \frac{1}{N^2}\sum_{i=1}^{N} \mathrm{Var}(f(X_i)) = \frac{\sigma^2}{N}, \quad \sigma^2 \equiv \mathrm{Var}_p(f(X)) $$
と書けます。標準誤差 (standard error) はこの平方根で $\sigma/\sqrt{N}$ となり、推定誤差は $1/\sqrt{N}$ の速度で減衰します。重要なのは、この式に次元 $d$ がまったく現れていないことです。1次元だろうが100次元だろうが、$\sigma^2$ が有限であれば誤差は同じ速さで減ります。これが冒頭で述べた「次元の呪いの回避」の正体です。
中心極限定理と信頼区間
実用上は「推定値の不確かさ」を区間で示したいことが多いので、中心極限定理 (CLT) を当てはめます。$\sigma^2 < \infty$ ならば $N \to \infty$ で
$$ \sqrt{N}\,(\hat{I}_N – I) \xrightarrow{d} \mathcal{N}(0, \sigma^2) $$
が成り立ちます。実際の計算では $\sigma^2$ も標本から推定するしかなく、
$$ \hat{\sigma}^2 = \frac{1}{N-1}\sum_{i=1}^{N}\big(f(X_i) – \hat{I}_N\big)^2 $$
を用いて 95% 信頼区間を
$$ \hat{I}_N \pm 1.96 \cdot \frac{\hat{\sigma}}{\sqrt{N}} $$
と構成します。この信頼区間が「ただ平均を出す」モンテカルロ法を、決定論的な数値積分と一線を画す道具にしています。決定論的アルゴリズムでは誤差の上界を出すのに高次の滑らかさが必要ですが、モンテカルロ法では分散一つで誤差を見積もれるのです。
次元の呪いを回避する理由
なぜ次元によらないのか、もう少し噛み砕いておきます。決定論的なグリッド法は「空間全体を覆う」必要があるため、点数が $K^d$ で増えます。一方、モンテカルロ法は分布 $p(x)$ が高密度なところに自然とサンプルが集中します。10次元正規分布の「中心 $1\sigma$ 球内」の体積は全体のごく一部ですが、$p$ からサンプルすればその大部分のサンプルがその中に落ちます。「測度の集中」を逆手に取って、関心のある領域だけを効率的に探索する——これがモンテカルロ法の本質です。
ただし、これには重大な但し書きがあります。$\sigma^2 = \mathrm{Var}_p(f)$ が有限でなければ収束が崩れ、有限であっても $\sigma^2$ が極端に大きければ実用的な精度に達するまでに天文学的な $N$ が必要になります。次節で扱う重点サンプリングは、まさにこの「$\sigma^2$ が大きすぎる問題」に対する処方箋です。
重点サンプリングの定式化
サンプリング元を取り替える
素朴MCでは「目的の分布 $p$ から直接サンプルできる」ことを暗黙に仮定していました。しかし現実には、$p$ から直接引けない（例: 正規化定数が未知のベイズ事後分布）、引けても効率が悪い（例: 関心領域が裾にある稀少事象）ことが多々あります。
ここで重点サンプリングの基本トリックは「別の分布 $q$ からサンプルして、重みで補正する」というものです。$q(x) > 0$ となる領域が $p(x)f(x) \neq 0$ の領域を含んでいる限り（$p \cdot f$ の台が $q$ の台に含まれる、いわゆる絶対連続性の条件）、次の恒等式が成り立ちます。
$$ I = \int f(x) p(x)\, dx = \int f(x) \frac{p(x)}{q(x)} q(x)\, dx = \mathbb{E}_q\!\left[f(X)\,w(X)\right] $$
ここで $w(x) \equiv p(x)/q(x)$ は重要度重み (importance weight) と呼ばれます。両辺がぴったり等しいのは、$dx$ について積分するときに $q(x)/q(x) = 1$ を挟むという単純な事実から来ています。新しい推定量は
$$ \hat{I}^{\mathrm{IS}}_N = \frac{1}{N}\sum_{i=1}^{N} f(X_i) w(X_i), \quad X_i \sim q $$
となります。これも不偏で、$\mathbb{E}_q[\hat{I}^{\mathrm{IS}}_N] = I$ が線形性からただちに従います。
分散はどう変わるか
不偏性は嬉しいですが、本当に効くのは分散です。重点サンプリング推定量の分散を計算してみましょう。
$$ \mathrm{Var}_q\!\left(\hat{I}^{\mathrm{IS}}_N\right) = \frac{1}{N}\,\mathrm{Var}_q\!\big(f(X)w(X)\big) = \frac{1}{N}\left( \mathbb{E}_q[f^2 w^2] – I^2 \right) $$
ここで第1項を $p$ のもとでの期待値に書き換えると、
$$ \mathbb{E}_q[f(X)^2 w(X)^2] = \int f(x)^2 \frac{p(x)^2}{q(x)^2} q(x)\, dx = \int f(x)^2 \frac{p(x)^2}{q(x)}\, dx $$
となり、最終的に
$$ \boxed{\;\mathrm{Var}_q\!\left(\hat{I}^{\mathrm{IS}}_N\right) = \frac{1}{N}\left( \int \frac{f(x)^2 p(x)^2}{q(x)}\, dx – I^2 \right)\;} $$
が得られます。この式から重要な示唆が読み取れます。$q(x)$ が分母にあるため、$f(x)p(x)$ が大きいのに $q(x)$ が小さい点があると、被積分関数が爆発し、分散が暴れます。逆に、$q(x)$ を $|f(x)|p(x)$ の形に近くなるように選べば、分子と分母が打ち消し合って分散はぐっと抑えられます。「提案分布の形をどうデザインするか」がIS成功の鍵だと、この式はずばり告げています。
自己正規化重点サンプリング (Self-Normalized IS, SNIS)
実用上、目的分布 $p$ や提案 $q$ の正規化定数が分からないことがよくあります。とくにベイズの事後分布 $p(\theta | D) = \tilde{p}(\theta)/Z$ では分母 $Z$（証拠）が高次元積分で、まさにこれを推定したくて重点サンプリングを使うのに、重み $w = p/q$ を計算するためには $Z$ が必要——という循環が生じます。
これを回避するのが自己正規化重点サンプリングです。未正規化の重み $\tilde{w}(x) = \tilde{p}(x)/\tilde{q}(x)$ を使い、
$$ \hat{I}^{\mathrm{SNIS}}_N = \frac{\sum_{i=1}^{N} f(X_i) \tilde{w}(X_i)}{\sum_{i=1}^{N} \tilde{w}(X_i)} $$
と定義します。分子分母にそれぞれ正規化定数比 $Z_q/Z_p$ がかかって相殺するため、未正規化のまま使えます。代償として推定量はバイアスを持ちますが、バイアスは $O(1/N)$ で消えていき、分散の $O(1/\sqrt{N})$ よりも速く消えるため、$N$ が十分大きければ無視できます。実装の手軽さから、現代のベイズ計算ではSNISが標準的に使われます。
提案分布の良し悪しが分散を支配するというのは分かりました。では、どんな $q$ が「最良」なのでしょうか。次のセクションで、理論上の最適形と、それを近似的に評価する指標を導きます。
最適提案分布とESS
分散をゼロにする提案分布
分散の式
$$ \mathrm{Var}_q\big(f(X)w(X)\big) = \mathbb{E}_q[f^2 w^2] – I^2 $$
を最小化する $q$ を考えます。$f(x) \geq 0$ と仮定して、ラグランジュ未定乗数法で「$\int q(x)\, dx = 1$ の制約下で $\mathbb{E}_q[f^2 w^2]$ を最小化する」問題を解くと、最適解は
$$ q^*(x) = \frac{f(x) p(x)}{\int f(x’) p(x’)\, dx’} = \frac{f(x) p(x)}{I} $$
となります。一般の $f$ では $|f|$ を使って
$$ \boxed{\;q^*(x) = \frac{|f(x)|\, p(x)}{\int |f(x’)|\, p(x’)\, dx’}\;} $$
が最適提案分布です。$f \geq 0$ のとき、この $q^*$ を実際に分散の式に代入すると、$f(x)w(x) = f(x)\cdot p(x)/q^*(x) = I$ となり、すべてのサンプルが同じ値 $I$ を返すため分散はぴったりゼロになります。1サンプルで真値が分かる、夢のような提案分布です。
ただし夢には罠があります。$q^*$ の正規化定数こそ求めたかった $I$ そのものであり、$q^*$ を構成できる時点で問題はすでに解けています。実用上は「$q^*$ にできるだけ似た形」の扱いやすい分布（多変量ガウス、$t$ 分布、混合正規など）を提案分布に取ることになります。この近似の品質をどう測るかが次の課題です。
青が目的分布 $p$、緑が被積分関数 $f(x)=e^{x/2}$、赤が最適提案 $q^*(x)\propto |f(x)|p(x)$ です。$f$ が大きい側（$x$ が正の方向）に $q^*$ が引っ張られているのが分かります。最適提案は「$f$ が大きくて $p$ も非ゼロな領域」にサンプルを集中させ、1サンプルあたりの情報量を最大化しています。
有効サンプル数 (Effective Sample Size, ESS)
提案分布の良し悪しを実際のサンプルから判定する指標が有効サンプル数 (ESS) です。Kishの公式と呼ばれる古典的な定義は
$$ \boxed{\;\mathrm{ESS} = \frac{\left(\sum_{i=1}^{N} w_i\right)^2}{\sum_{i=1}^{N} w_i^2}\;} $$
で、$1 \leq \mathrm{ESS} \leq N$ の範囲を取ります。極端なケースを見ましょう。
- 全ての重みが等しい（$w_1 = \dots = w_N$）とき、$\mathrm{ESS} = N$。これは「提案分布が目的分布と一致しているので、サンプル全てが有効」状態です。
- 1つの重みだけ大きく他がゼロのとき、$\mathrm{ESS} = 1$。これは「実質1サンプルしか効いていない」破滅的状態です。
ESSの意味は「同じ精度の素朴MC推定量を作るのに何個の独立サンプルが必要か」という換算量です。$N=10000$ で $\mathrm{ESS}=50$ なら、200倍のサンプルを引いたのに実質50サンプル分の情報しか取れていない、ということ。重点サンプリングを使う際は、必ずESSをモニターし、$\mathrm{ESS}/N$ が0.1未満になったら提案分布の再設計を考えるのが定石です。
左（良い提案、$q \approx p$）では重みがほぼ均一に分布し、ESS/N ≈ 0.8 を保っています。右（悪い提案、$q$ が狭い）では重みが一点に偏り、ESS/N が急落します。この「重みの均一性」がESSの本質で、ESS/Nを見るだけで提案分布の適切さを即座に診断できます。実装では推定値と同時にESSも常に出力する習慣を付けましょう。
重み崩壊 (Weight Degeneracy)
ESSが極端に小さくなる現象を重み崩壊 (weight degeneracy) と呼びます。これは次元が高くなるほど顕著になります。直感的には、$d$ 次元空間で2つの分布の「重なり」が指数的に小さくなるためです。
具体例として、目的が $d$ 次元標準正規 $\mathcal{N}(0, I_d)$、提案が分散の少し違う $\mathcal{N}(0, \sigma^2 I_d)$ の場合、重みの分散は
$$ \mathrm{Var}_q(w) \propto \left(\frac{\sigma^2}{2\sigma^2 – 1}\right)^{d/2} – 1 $$
のように $d$ について指数的に増えていきます（$\sigma^2 > 1/2$ のとき）。$d=100$ ともなれば、$\sigma^2 = 1.1$ という1割の差でもESSが急激に劣化します。これが「素朴ISは高次元で機能しない」と言われる所以で、後述するAnnealed Importance SamplingやSequential Monte Carloといった発展手法が必要になる理由でもあります。
バイアス・バリアンス分解
実用上、提案分布のチューニングはバイアスとバリアンスのトレードオフになります。SNISの平均二乗誤差 (MSE) は
$$ \mathrm{MSE}(\hat{I}^{\mathrm{SNIS}}_N) = \mathrm{Var}(\hat{I}^{\mathrm{SNIS}}_N) + \mathrm{Bias}(\hat{I}^{\mathrm{SNIS}}_N)^2 \approx \frac{C_1}{N} + \frac{C_2}{N^2} $$
のように分解できます。$N \to \infty$ では分散項が支配的なので、長期的には分散低減が最優先です。短いランで提案分布を試行錯誤するときはバイアスが顔を出すことに注意します。
最適提案分布の理論とESSによる診断が揃いました。次は、重点サンプリングが爆発的に効く具体例——稀少事象推定——を見ていきます。
レアイベント推定
なぜ素朴MCはレアイベントに弱いのか
レアイベント (rare event) とは、確率 $p^* \ll 1$ で起こる事象 $A$ のことです。$P(X \in A)$ を素朴MCで推定すると、
$$ \hat{p}^*_N = \frac{1}{N}\sum_{i=1}^{N} \mathbb{1}[X_i \in A], \quad X_i \sim p $$
となります。$\mathbb{1}[\cdot]$ は指示関数で、事象 $A$ が起きたサンプルだけが1にカウントされます。この推定量の相対誤差は
$$ \frac{\sqrt{\mathrm{Var}(\hat{p}^*_N)}}{p^*} = \sqrt{\frac{p^*(1-p^*)}{N (p^*)^2}} \approx \frac{1}{\sqrt{N p^*}} $$
となり、$p^* = 10^{-6}$ で 1% の相対誤差を出すには $N \approx 10^{10}$ サンプル必要です。これは現実的とは言えません。
問題の本質は「$A$ に落ちるサンプルが極めて少ない」ことです。たとえば $X \sim \mathcal{N}(0,1)$ で $P(X > 4)$ を計算しようとすると、$4\sigma$ より外のサンプルは平均すると $\sim 3 \times 10^{-5}$ の確率でしか現れず、$N=10^6$ サンプルでも数十個しか拾えません。
IS による打開策
重点サンプリングなら、$A$ に多くのサンプルを意図的に落とすような $q$ を選んで、その分を重みで補正できます。$P(X > 4)$ の例なら、$X$ を平均4の正規分布 $\mathcal{N}(4, 1)$ から引けば、ほぼ半分のサンプルが $A = \{X > 4\}$ に落ちます。重みは
$$ w(x) = \frac{p(x)}{q(x)} = \frac{(2\pi)^{-1/2}\exp(-x^2/2)}{(2\pi)^{-1/2}\exp(-(x-4)^2/2)} = \exp\!\left(-\frac{x^2 – (x-4)^2}{2}\right) = \exp(-4x + 8) $$
となります。$x = 4$ 付近では $w \approx e^{-8} \approx 3.4 \times 10^{-4}$ という小さな値で、これによって「$\mathcal{N}(0,1)$ では $x=4$ 周辺はほとんど起こらない」事実が正しく重みに反映されます。シミュレーションすると、$N=10^4$ サンプルで素朴MCの $N=10^7$ サンプルと同等以上の精度が得られます。1000倍以上の効率改善です。
提案分布の選び方の経験則
レアイベント推定で提案分布を選ぶときの実用的な指針をまとめておきます。
- 指数傾斜 (Exponential tilting): 目的が指数族 $p(x) \propto e^{-\phi(x)}$ なら、$q(x) \propto e^{-\phi(x) + \lambda^\top x}$ のように線形項で傾斜させると、平均が $A$ の方向にシフトします。クラメール定理の大偏差原理から自然に出てくる選び方です。
- 混合分布: 目的の裾と提案の裾を両方カバーするよう、$q = (1-\alpha) p + \alpha q_{\mathrm{tail}}$ のように混ぜます。重み爆発を防げます。
- 適応的IS (Adaptive IS): 試行錯誤で得たサンプルから $q$ のパラメータを更新していく方法。Cross-Entropy法、Population Monte Carloなどが代表例。
応用例: 金融VaR、稀少事象の確率推定
レアイベントISは産業界で広く使われています。
- 金融VaR (Value at Risk): 「翌日損失が $L$ を超える確率」を見積もる際、$L$ が大きいほどレアです。指数傾斜やヘビーテイル提案分布を用いたISが用いられます。
- 移動体のニアミス確率: 2機の航空機や車両の位置予測の不確かさをガウス分布で表し、接近事象（最接近距離 < 安全マージン）の確率を計算します。素朴MCでは接近が稀すぎて推定が荒れるため、最接近点周辺を厚くサンプリングする提案分布が使われます。
- 半導体故障率推定: VLSIチップの動作不良率は $10^{-9}$ オーダー。素朴MCでは絶望的なので、製造ばらつき空間に偏らせたISが標準的に使われます。
理論はここまで揃いました。次のセクションで、これらの主張をPythonで一つひとつ実証していきましょう。
Python実装 — MC vs IS の効率比較
1次元積分での分散比較
まず素朴MCとISの基本的な動きを、解析解が分かる積分で比較します。$\int_{-\infty}^{\infty} x^2 \phi(x; 0, 1)\, dx = 1$ という単純な期待値（標準正規分布の2次モーメント）を例に取ります。
import numpy as np
np.random.seed(0)
N = 10000
# 素朴MC: 標準正規からサンプル
x_p = np.random.randn(N)
est_crude = np.mean(x_p ** 2)
var_crude = np.var(x_p ** 2) / N
# IS: 提案分布として N(0, 2^2) を使う（裾を厚くする）
sigma_q = 2.0
x_q = sigma_q * np.random.randn(N)
log_p = -0.5 * x_q ** 2 - 0.5 * np.log(2 * np.pi)
log_q = -0.5 * (x_q / sigma_q) ** 2 - 0.5 * np.log(2 * np.pi * sigma_q ** 2)
w = np.exp(log_p - log_q)
est_is = np.mean((x_q ** 2) * w)
var_is = np.var((x_q ** 2) * w) / N
print(f"真値 : 1.0")
print(f"素朴MC : {est_crude:.4f} 分散 = {var_crude:.4e}")
print(f"IS(σ_q=2) : {est_is:.4f} 分散 = {var_is:.4e}")
print(f"分散比 IS/MC: {var_is/var_crude:.3f}")
実行すると、IS推定量の分散が素朴MCの約2.5倍になることが分かります。「$x^2$ の期待値」を計算するだけなら素朴MCで十分で、わざわざ裾を広げると重みのばらつきが分散を悪化させるのです。これは重要な教訓で——何でもかんでもISを使えばよいわけではない——提案分布が悪いと素朴MCより悪化します。次の節でこの逆、ISが圧勝するケースを見ます。
2000回の試行で得た推定値のヒストグラムです。素朴MCの分布は真値を中心に広く散らばり、ゼロ付近に多くの試行が集中します（サンプルがヒットしないため）。一方ISの分布は真値のほぼ直下に鋭く集中しており、標準偏差の比が約10〜100倍になることが視覚的に確認できます。
レアイベント $P(X>4)$ の1000倍効率化
ISが本領を発揮するのはレアイベントです。$P(X > 4)$ を素朴MCとISで比較します。
import numpy as np
from scipy.stats import norm
np.random.seed(1)
N = 100000
threshold = 4.0
true_prob = 1 - norm.cdf(threshold)
print(f"真値 P(X > {threshold}) = {true_prob:.4e}")
# 素朴MC
x_p = np.random.randn(N)
hits = (x_p > threshold).astype(float)
est_crude = hits.mean()
se_crude = hits.std(ddof=1) / np.sqrt(N)
print(f"素朴MC : 推定値={est_crude:.4e}, SE={se_crude:.2e}, 衝突回数={int(hits.sum())}")
# IS: N(4, 1) を提案分布に使う
shift = 4.0
x_q = shift + np.random.randn(N)
w = np.exp(-shift * x_q + 0.5 * shift ** 2) # p(x)/q(x) = exp(-shift*x + shift^2/2)
contrib = (x_q > threshold).astype(float) * w
est_is = contrib.mean()
se_is = contrib.std(ddof=1) / np.sqrt(N)
print(f"IS(平均4) : 推定値={est_is:.4e}, SE={se_is:.2e}")
print(f"標準誤差比 SE_MC/SE_IS = {se_crude/se_is:.1f} 倍")
print(f"実質的なサンプル効率比 = {(se_crude/se_is)**2:.0f} 倍")
このコードを実行すると、素朴MCでは $N=10^5$ サンプルのうち平均すると数個しか $x > 4$ に該当しません。たまたま0回になれば推定値もゼロです。一方ISでは半数近くが該当領域に落ち、推定値が安定します。標準誤差比は数十倍、効率比（その2乗）は1000倍前後になります。同じ精度を素朴MCで得るには、サンプル数を1000倍にする必要があるという意味です。理論で予測した「レアイベントでISが劇的に効く」が数値で確認できました。
左図（素朴MC）では閾値 $x=4$ を超えるサンプルがほとんど現れず、推定精度が極めて低いことが分かります。右図（IS）では提案分布を $N(4,1)$ にシフトすることで、約半数のサンプルが閾値を超え、確率の推定に直接貢献します。サンプル数は同じでも、稀少域への集中度が劇的に異なる点が視覚的に明確です。
多次元での分散と次元依存性
ISの効果は次元によって大きく変わります。$d$ 次元標準正規からの $\mathbb{E}[\|X\|^2] = d$ をターゲットに、提案分布を $\mathcal{N}(0, 1.5^2 I_d)$ として、$d$ ごとに重みの分散とESSを見てみます。
import numpy as np
def is_estimate(d, N=10000, sigma_q=1.5, seed=0):
"""d次元正規での重みESS計算"""
rng = np.random.default_rng(seed)
x = sigma_q * rng.standard_normal((N, d))
# log p(x) - log q(x) を成分ごとに計算（数値安定のためlogで処理）
log_p = -0.5 * np.sum(x ** 2, axis=1) - 0.5 * d * np.log(2 * np.pi)
log_q = -0.5 * np.sum((x / sigma_q) ** 2, axis=1) - 0.5 * d * np.log(2 * np.pi * sigma_q ** 2)
log_w = log_p - log_q
# 数値安定化
log_w -= log_w.max()
w = np.exp(log_w)
ess = (w.sum() ** 2) / (w ** 2).sum()
f = np.sum(x ** 2, axis=1)
est = (f * w).sum() / w.sum()
return est, ess
for d in [1, 5, 10, 20, 50, 100]:
est, ess = is_estimate(d, N=10000, sigma_q=1.5)
print(f"d={d:3d}: 推定値={est:8.2f} (真値={d:3d}) ESS={ess:8.1f} / 10000 ESS率={ess/10000:.4f}")
このコードを実行すると、低次元（$d=1$）ではESSが数千あるのに対し、$d=50$ や $d=100$ ではESSが1桁台にまで激減します。これがまさに重み崩壊で、提案分布のわずかな選択ミス（$\sigma_q = 1.5$ という穏当な値ですら）が高次元では致命的になります。逆に言えば、「ISが効くかどうか」は次元と提案分布の質に強く依存し、無策に高次元へ適用すると失敗するということです。実用ではAnnealed Importance Sampling (AIS) やSequential Monte Carloで段階的に分布を変形する戦略が取られます。
各曲線が提案分布の標準偏差 $\sigma_q$ を表しています。$\sigma_q=1.1$（目的とほぼ同じ）でさえ $d=50$ 付近でESS/Nが0.1を大幅に下回ります。$\sigma_q$ が大きくなるほど崩壊が加速し、$d=100$ では全ての提案でESSがほぼゼロになることが分かります。縦軸が対数スケールである点に注目——次元増加によるESSの低下は本当に指数的です。
ESSの劣化過程を可視化
提案分布の「ずれ具合」を変えながらESSがどう劣化するか、グラフで見てみましょう。
import numpy as np
import matplotlib.pyplot as plt
def ess_curve(d, sigma_qs, N=10000, seed=42):
rng = np.random.default_rng(seed)
ess_list = []
for sq in sigma_qs:
x = sq * rng.standard_normal((N, d))
log_p = -0.5 * np.sum(x ** 2, axis=1)
log_q = -0.5 * np.sum((x / sq) ** 2, axis=1) - d * np.log(sq)
lw = log_p - log_q
lw -= lw.max()
w = np.exp(lw)
ess = (w.sum() ** 2) / (w ** 2).sum()
ess_list.append(ess / N)
return ess_list
sigma_qs = np.linspace(0.7, 2.0, 30)
plt.figure(figsize=(9, 5))
for d in [1, 5, 10, 20, 50]:
plt.plot(sigma_qs, ess_curve(d, sigma_qs), 'o-', label=f'd={d}', markersize=4)
plt.axvline(1.0, color='k', ls='--', alpha=0.5, label='σ_q = σ_p')
plt.xlabel('Proposal scale σ_q (target σ_p = 1)')
plt.ylabel('ESS / N')
plt.yscale('log')
plt.title('ESS degradation vs proposal mismatch')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ess_degradation.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、いくつかの重要な傾向が読み取れます。第一に、$\sigma_q = 1$ では全次元で ESS/N = 1（提案分布が目的分布と一致するから当然）。第二に、$\sigma_q$ が1から離れるにつれESSが急速に下がり、その下がり方が次元 $d$ について指数的に強まります。第三に、提案分布が目的より狭い（$\sigma_q < 1$）方が、広い（$\sigma_q > 1$）よりも遥かに危険です——これは $w \propto e^{x^2(\sigma_q^{-2}-1)/2}$ が裾で爆発するためで、「目的より広く取る」が安全側のヒューリスティクスになっています。実装でISを使うときは、まずESSをこのようにスキャンして提案分布のロバスト性を確認するのが定石です。
左図（分散）では $\sigma_q=1$ 付近でほぼ最小となり、そこから離れると急速に分散が増大します。特に $\sigma_q < 1$ 側の増大が急峻で、重みが爆発して分散が無限大に近づく様子が対数スケールで確認できます。右図（ESS）は左の鏡像で、$\sigma_q=1$ でESS最大、狭い側で急落します。「少し広め」が最も安全な選択であることをこの2図が定量的に示しています。
信頼区間の幅を比較
最後に、$P(X>4)$ 推定での95%信頼区間を素朴MCとISでプロットし、誤差の縮小を視覚的に確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
true_p = 1 - norm.cdf(4.0)
Ns = np.logspace(2, 6, 20).astype(int)
def crude_estimate(N, seed):
rng = np.random.default_rng(seed)
x = rng.standard_normal(N)
h = (x > 4.0).astype(float)
return h.mean(), h.std(ddof=1) / np.sqrt(N) if N > 1 else np.nan
def is_estimate_rare(N, seed):
rng = np.random.default_rng(seed)
x = 4.0 + rng.standard_normal(N)
w = np.exp(-4.0 * x + 8.0)
c = (x > 4.0).astype(float) * w
return c.mean(), c.std(ddof=1) / np.sqrt(N) if N > 1 else np.nan
est_c, se_c = [], []
est_i, se_i = [], []
for N in Ns:
a, b = crude_estimate(N, seed=2025)
est_c.append(a); se_c.append(b)
a, b = is_estimate_rare(N, seed=2025)
est_i.append(a); se_i.append(b)
est_c = np.array(est_c); se_c = np.array(se_c)
est_i = np.array(est_i); se_i = np.array(se_i)
plt.figure(figsize=(10, 5))
plt.fill_between(Ns, est_c - 1.96 * se_c, est_c + 1.96 * se_c, alpha=0.3, label='Crude MC 95% CI')
plt.plot(Ns, est_c, 'o-', label='Crude MC mean')
plt.fill_between(Ns, est_i - 1.96 * se_i, est_i + 1.96 * se_i, alpha=0.3, label='IS 95% CI')
plt.plot(Ns, est_i, 's-', label='IS mean')
plt.axhline(true_p, color='k', ls='--', label=f'True = {true_p:.2e}')
plt.xscale('log'); plt.yscale('symlog', linthresh=1e-6)
plt.xlabel('Sample size N')
plt.ylabel('Estimated P(X > 4)')
plt.title('Confidence intervals: Crude MC vs Importance Sampling')
plt.legend(loc='upper right', fontsize=9)
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('ci_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、3つの本質的な事実が読み取れます。第一に、素朴MCの信頼区間は $N$ が $10^4$ 程度までゼロを含み、推定値自体も真値からかなり外れます（衝突回数が0または1のときの統計は壊れます）。第二に、IS推定量の信頼区間は $N \sim 10^2$ 程度でも真値を捉え、$N$ とともに $1/\sqrt{N}$ で順調に縮みます。第三に、両者の信頼区間幅を同じ $N$ で比べると、IS の方が桁違いに狭く、レアイベントでのISの威力が一目瞭然です。理論で導いた分散低減が、信頼区間という実用的な形で実証できました。
縦軸は対数スケールで、赤い破線が真値 $P(X>4)\approx 3.2\times 10^{-5}$ です。青の素朴MCは $N=10^4$ でもまだ信頼区間が真値を大きく外れ、$N=10^5$ でやっと安定し始めます。緑のISは $N=100$ 程度から既に真値付近を捉え、信頼区間の幅も同じ $N$ で素朴MCの100分の1程度です。両者の収束速度の差が対数スケールでも歴然としており、レアイベント推定でISを使うことの必然性を視覚化しています。
Self-normalized IS でベイズ事後の期待値を求める
実用上は正規化定数が未知の事後分布が頻出します。SNISでこれを扱うミニ実装を見ましょう。ターゲットを「未正規化事後」 $\tilde{p}(\theta) = \exp(-\theta^2/2) \cdot (1 + 0.5\sin(3\theta))$（標準正規 × 周期的な揺らぎ）と置き、$\mathbb{E}_p[\theta^2]$ を求めます。
import numpy as np
import matplotlib.pyplot as plt
def log_tilde_p(theta):
"""未正規化対数事後（正規化定数は不要）"""
return -0.5 * theta ** 2 + np.log(1.0 + 0.5 * np.sin(3.0 * theta))
def log_q(theta, sigma=1.5):
"""提案分布: 平均0・標準偏差sigma のガウス"""
return -0.5 * (theta / sigma) ** 2 - 0.5 * np.log(2 * np.pi * sigma ** 2)
rng = np.random.default_rng(7)
N = 50000
sigma_q = 1.5
theta = sigma_q * rng.standard_normal(N)
# 未正規化重み（数値安定化のためlog空間で）
log_w = log_tilde_p(theta) - log_q(theta, sigma_q)
log_w -= log_w.max()
w = np.exp(log_w)
# Self-normalized IS 推定量
f = theta ** 2
est_snis = np.sum(f * w) / np.sum(w)
ess = (w.sum() ** 2) / (w ** 2).sum()
Z_est = w.mean() * np.exp(log_w.max() + 0.5 * np.log(2 * np.pi * sigma_q ** 2))
# 参考: 数値積分でのground truth
from scipy.integrate import quad
Z_true, _ = quad(lambda t: np.exp(log_tilde_p(t)), -10, 10)
num, _ = quad(lambda t: t ** 2 * np.exp(log_tilde_p(t)), -10, 10)
true_val = num / Z_true
print(f"SNIS推定 E[θ²] = {est_snis:.4f}")
print(f"数値積分の真値 = {true_val:.4f}")
print(f"ESS / N = {ess/N:.3f}")
print(f"正規化定数 Z 推定 = {Z_est:.4f} 真値 = {Z_true:.4f}")
このコードを実行すると、SNIS推定値と数値積分の真値はほぼ一致し、相対誤差は1%程度で収まります。ESS/Nは0.5〜0.8 程度で、提案分布が事後を十分カバーできていることを示します。正規化定数 $Z$ も合わせて推定できる点が重要で、これがベイズモデル選択（ベイズ因子の計算）に直結します。「未正規化分布のままで使える」というSNISの実用上の利点が、この例で明確に確認できました。
左図では未正規化事後 $\tilde{p}(\theta)$（青）と提案分布 $q=N(0, 1.5^2)$（橙の破線）を重ねています。提案が目的を包んでいることが確認でき、ESS/N ≈ 0.74 という高い値もこの良好な重なりを反映しています。右図はSNIS重みのヒストグラムで、ほぼゼロ付近に集中し少数の大きな重みが存在する様子が見えます。この非均一性がバイアスの源ですが、$N=20000$ では実用上問題ない精度（誤差0.3%程度）が得られています。
応用とMCMCへの接続
ベイズ周辺尤度の計算
モデル比較で重要な周辺尤度（証拠）
$$ Z = p(D) = \int p(D | \theta) p(\theta)\, d\theta $$
は典型的な高次元積分です。提案分布 $q(\theta)$ を用いて
$$ Z = \mathbb{E}_q\!\left[\frac{p(D | \theta)p(\theta)}{q(\theta)}\right] $$
と評価できます。しかし、事前 $p(\theta)$ と事後 $p(\theta|D) \propto p(D|\theta)p(\theta)$ の重なりが小さいと（高次元では典型的）、$q = p$ では重み崩壊が起きます。$q$ をラプラス近似や変分近似で事後の形に合わせる工夫が要ります。ブリッジサンプリングやパスサンプリングは、$Z$ の計算に特化したIS変種で、ベイズモデル選択の現代的標準です。
強化学習のオフポリシー評価
強化学習で「収集したデータは振る舞いポリシー $\mu$ で取られたが、評価したいのは目標ポリシー $\pi$ である」状況がよくあります。期待リターン
$$ J(\pi) = \mathbb{E}_\pi[R] = \mathbb{E}_\mu\!\left[\frac{\pi(a|s)}{\mu(a|s)} R\right] $$
がIS推定量です。1ステップあたりの重要度比 $\rho_t = \pi(a_t|s_t)/\mu(a_t|s_t)$ をエピソード全体で積算するため、エピソード長が長いと重み崩壊が起きやすく、Per-Decision IS、Weighted IS、Doubly Robust推定などの工夫が研究されています。
変分推論との関係
変分推論 (VI) では、扱いやすい分布族 $q_\phi$ で事後を近似し、ELBO
$$ \mathcal{L}(\phi) = \mathbb{E}_{q_\phi}\!\left[\log\frac{p(\theta, D)}{q_\phi(\theta)}\right] $$
を最大化します。最適な $q_\phi$ は「事後への KL ダイバージェンス最小」を満たし、これがIS提案分布の良い候補にもなります。実際、Importance Weighted Autoencoder (IWAE) はVIの $q_\phi$ をIS提案分布として使い、複数サンプルの重みでELBOを多重サンプリングする方法で、変分近似の偏りを定量的に補正します。
MCMCへの自然な接続
提案分布 $q$ を「固定された分布」ではなく「現在の状態 $\theta_t$ に依存する条件付き分布 $q(\theta’|\theta_t)$」に拡張すると、Metropolis-Hastings法に到達します。新しい候補 $\theta’$ を $q(\theta’|\theta_t)$ から提案し、
$$ \alpha = \min\!\left(1, \frac{p(\theta’)q(\theta_t|\theta’)}{p(\theta_t)q(\theta’|\theta_t)}\right) $$
の確率で受理する——これは「IS重みに基づく受理確率」と解釈できます。ISが「独立な提案」を扱うのに対し、MCMCは「逐次的に提案を更新する」枠組みで、両者は提案分布の使い方こそ違うものの、根本では同じ重み補正の発想に立脚しています。
左図のMH法では二峰性の目標分布（赤）を逐次提案で再現できています（受理率 ≈ 0.55）。右図のISでは幅広い正規提案から重み付きで同じ分布を再現しており、ESS/N ≈ 0.36 と比較的健全です。1次元の二峰性程度であればどちらも機能しますが、高次元になるとISのESSが急落するため、MCMCの方が実用的です。ISとMCMCは「重みで分布を修正する」というコアアイデアを共有しており、用途に応じて使い分けることが重要です。
Annealed Importance Sampling (AIS)
高次元でISが崩壊する問題への決定的な解決策がAnnealed ISです。事前 $p_0$ と目的事後 $p_K$ の間に中間分布 $p_k \propto p_0^{1-\beta_k} p_K^{\beta_k}$（$0 = \beta_0 < \dots < \beta_K = 1$）を挟み、各段でMCMCを少しずつ実行しながら重みを蓄積します。物理学のシミュレーテッドアニーリングと類似の構造で、温度を冷ましながら徐々に目的分布へ近づけます。AISは変分推論より遥かに高品質な周辺尤度推定を与え、深層生成モデルの評価でデファクト標準になっています。
まとめ
本記事では、高次元期待値計算の中核技術であるモンテカルロ法と重点サンプリングを、理論の導出から実装・応用まで通して解説しました。
- 素朴モンテカルロ法: 積分を期待値 $I = \mathbb{E}_p[f]$ に書き換え、$\hat{I}_N = \frac{1}{N}\sum f(X_i)$ で推定する。不偏で、分散は $\sigma^2/N$、誤差は次元に依らず $O(1/\sqrt{N})$。CLTにより信頼区間 $\hat{I}_N \pm 1.96 \hat{\sigma}/\sqrt{N}$ が引ける。
- 重点サンプリング: 提案分布 $q$ からサンプルして $f(x)p(x)/q(x)$ を平均する。分散 $\mathrm{Var}_q(\hat{I}^{\mathrm{IS}}) = \frac{1}{N}(\int f^2 p^2 / q\,dx – I^2)$ は $q$ の選び方で激変する。
- 最適提案分布とESS: $q^*(x) \propto |f(x)|p(x)$ で分散がゼロ（実用不能）。実装ではESS $= (\sum w)^2/\sum w^2$ で提案分布の品質を診断する。
- レアイベント推定: 素朴MCの相対誤差 $1/\sqrt{Np^*}$ は $p^* \ll 1$ で破綻するが、ISは正規分布の例で1000倍の効率改善を達成する。
- 罠と対策: 重み崩壊（高次元で深刻）、self-normalized IS（正規化定数未知でも使える）、Annealed IS（高次元事後への決定打）。
- 応用と接続: ベイズ周辺尤度、強化学習オフポリシー評価、変分推論ELBO、そしてMetropolis-Hastings法への自然な拡張までを貫く一本の糸。
モンテカルロ法と重点サンプリングは、現代のベイズ計算・機械学習・計算物理学・金融工学・信頼性工学すべての基盤技術です。「サンプリングは積分より易しい」「重みで分布を取り替える」というシンプルな原理が、これほど広い領域を支えていることに驚かされます。
次のステップとして、以下の記事も参考にしてください。