シミュレーション：モンテカルロ法は遅いのではなく次元に鈍感【第34回】

はじめに

第32章はシミュレーションです。そして連載本編の最終回になります。

この連載では最初の回から、ほとんど毎回シミュレーションを使ってきました。「理論値と実測値を並べて確かめる」というやり方です。ところがその中身は、毎回こう書いていただけでした。

rng = np.random.default_rng(42)
x = rng.exponential(1.0, 100000)

rng.exponential の中で何が起きているのかを、私は一度も考えたことがありませんでした。 一様乱数しかない状態から、どうやって指数分布や正規分布が出てくるのか。第32章はその中身です。連載でずっと使ってきた道具の裏側を、最後に回収する回になります。

学び始める前の私は、この章を「テクニックの詰め合わせ」だと思っていました。逆関数法、棄却法、ボックス=ミュラー法、ブートストラップ法、並べ替え検定、分散減少法。名前が10個くらい並んでいて、それぞれ別に覚えるものだと。

ところが手を動かしてみて、いちばん効いたのは個々の手法ではなく章の構造でした。この章に出てくる技術は2種類しかありません。乱数を「作る」技術と、乱数を「使う」技術です。この2つを分けた瞬間に、10個の名前がきれいに2列に並びました。

そしてもう1つ、思い込みが崩れました。モンテカルロ法は速くも正確でもありません。 1変数の期待値を求めるだけなら、等間隔グリッドの数値積分に1億倍の精度差で負けます。実測しました。それでも使われるのは、精度で勝つからではなく「他の方法が全滅する場所でも生き残る」からでした。よく言われる「精度は $1/\sqrt{n}$ でしか上がらない」は、遅いという欠点ではなく「次元に鈍感」という長所だったのです。

この回で扱う用語

この章の地図：乱数を「作る」技術と「使う」技術

先に地図を出します。この章の理解はここで8割決まりました。

①の出力が②の入力になります。 そして「モンテカルロ法」という言葉は、本来②だけを指します。①は「乱数生成法」または「サンプリング法」です。

私はこの区別が付いておらず、棄却法とモンテカルロ法を同じものだと思っていました。混乱するのは自然で、というのも実際にまったく同じデータから両方が作れるからです。

正方形にダーツを400万本投げて、円の中に入ったかどうかを記録します。

• 緑の点の「個数」だけを使う → 円の面積比 $\pi/4$ の推定値が出る → $\pi = 3.142276$（真値 3.141593）。これがモンテカルロ法です。座標はもう捨ててよい
• 緑の点の「座標」を使う → 円内の一様分布に従う標本が手に入る。これが棄却法です。個数はどうでもよい

投げたダーツは1組しかありません。個数を使うか座標を使うかだけが違います。 そして欲しいものが「数値」か「標本」かで、呼び名が変わる。

以下、第1部で①を、第2部で②を扱います。

第1部　乱数を「作る」技術

逆関数法：累積分布関数を横向きに読む

手元にある道具は rng.random() ひとつ、つまり 0 と 1 のあいだを平らに埋める一様乱数だけ、という状況から出発します。ここから指数分布の乱数を作りたい。

累積分布関数 $F(x)$ は、いつもは横軸から縦軸に読みます。「$x=1$ のとき、それ以下になる確率は 0.632」。逆関数法がやるのは、この読み方をひっくり返すことだけです。縦軸に一様乱数 $u$ を置いて、曲線にぶつけて、横軸に落とす。

縦軸には 0.1 刻みで等間隔に $u$ を置いています。一様乱数なので、どの高さも等確率で出ます。ところが曲線を経由して横軸に落とすと、着地点は左に密集し右はまばらになります。 9本の矢印の間隔が、左では狭く右では広い。

この偏りが、そのまま指数分布のかたちです。

なぜ偏るのか：傾きが密度だから

左側で矢印が密集しているのは、そこで曲線の傾きが急だからです。傾きが急な区間では、狭い $x$ 幅のあいだに縦軸の広い範囲が対応します。縦軸の範囲が広い、つまりそこに落ちる $u$ が多い、つまりその $x$ の近くに乱数が多く生まれる。

そして累積分布関数の傾きは、微分すれば確率密度関数そのものです。

$F(x)$ の $x=0.2$ での傾きは 0.819、$x=2.5$ では 0.082。約10倍違います。右のパネルの $f(x)=e^{-x}$ の値を見ると、まさに同じ 0.819 と 0.082。「傾き」と「密度」は別の話ではなく、同じ数字の別名でした。

ここが原理のすべてです。「傾きが急なところに多く落ちる」という図の見た目が、そのまま「密度が高いところに多く落ちる」という確率の主張になっている。逆関数法は、密度の情報を累積分布関数の傾きという形で借りてきて、それを一様乱数の落下先の混み具合に翻訳する装置でした。

手続きは2行

指数分布なら $F(x) = 1 - e^{-x}$ を $u$ について解いて、$X = -\log(1-u)$ です。実際にやってみます。

乱数生成器の側は、指数分布のことを何も知りません。 rng.random() を呼んで $-\log(1-u)$ に通しただけです。

なぜ正しいのか：区間が1対1に対応している

「傾きが密度」という説明は直感的ですが、もう一段だけ厳密にしておきます。ここが分かると逆関数法は完全に手の内に入ります。

準備は2つだけです。

準備1。累積分布関数の縦軸は確率です。 だから縦軸上のある区間の「長さ」は、そのまま確率の値になっています。累積分布関数は「確率を長さに変換する装置」だと言えます。

準備2。一様乱数は「長さがそのまま確率になる」乱数です。 区間 $[a, b]$ に入る確率は $b - a$。定義から。

$x$ 軸側の帯は幅がバラバラ（0.5 / 1.0 / 1.5）ですが、$u$ 軸側の帯の長さはそれぞれの区間の確率そのものになっています。

あとは3段です。

1. $u$ が帯 $[F(1), F(2)]$ に落ちる確率は、その長さ = 0.2325（一様乱数の定義から）
2. 帯に落ちたら、$F$ は増加関数なので $F^{-1}(u)$ は $[1, 2]$ に入る。この例（連続で狭義に増加）では1対1に対応しているので、取りこぼしも重複もない
3. よって $X$ が $[1,2]$ に入る確率も 0.2325。そしてこれが、もともと欲しかった分布での $P(1 \le X \le 2)$ と同じ数字

式にすると1行です。

2番目の等号が「$F$ が単調非減少なので不等号の向きが保たれる」（厳密には次の節で述べる一般化逆関数の性質 $F^{-1}(u) \le x \iff u \le F(x)$）、3番目が「一様乱数が $[0, F(x)]$ に入る確率は長さ $F(x)$」。そして最後の $F(x)$ が、まさに目標の分布の累積分布関数です。到達しました。

実測で確認しました。 一様乱数200万個について「$u$ が帯に入ったか」と「$X$ が区間に入ったか」を1個ずつ比べたところ、200万個すべてで判定が一致（不一致 0 件）でした。

左2列が小数第5位まで完全一致しているのがポイントです。確率が近いのではなく、同じ事象を2つの軸から見ているだけでした。

だから「どんな分布でも」作れる

いま使った性質は「$F$ が単調非減少で右連続である」ことだけです。だから分布ごとに新しい理屈は要りません。

ただし1つだけ条件を正確にしておきます。累積分布関数が必ず満たすのは「単調非減少」と「右連続」で、狭義の単調増加ではありません。 離散分布や、密度が0の区間がある分布では、$F$ に水平な部分があって1対1になりません。そのときは逆関数を

と定義し直します（一般化逆関数）。こう定義すると $F^{-1}(u) \le x \iff u \le F(x)$ が一般に成り立つので、上の1行の証明はそのまま通ります。 狭義増加を仮定しなくてよいのが、この定義の効き目です。

実演します。同じ10個の一様乱数（0.07, 0.15, 0.28, … 0.95）を、4つの違う $F^{-1}$ に通しました。

4枚とも、灰色の水平線（$u$ の高さ）は同一です。違うのは曲線の形だけ。

読み取れることが3つあります。

• 曲線が立っている（密度が高い）ところでは着地点が密集し、寝ているところでは散らばる
• 離散分布（右下）は、累積分布関数が確率 $p_k$ の分だけ垂直に跳ぶので、その跳びの高さに落ちた $u$（幅 $p_k$ の区間）がすべて同じ値に写ることで実現される。連続分布と同じ手続きで離散分布も作れる
• 二山の混合分布では、$u = 0.28$ で $-1.250$ だったのが $u = 0.36$ で $+1.632$ に飛んでいます。$u = 0.30$ 付近が2つの山の谷間（曲線がほぼ水平）だからです。$u$ は等間隔なのに飛び地ができる

4つとも理論値と照合しました（各200万個）。

歪んだサイコロの各面の割合は 0.3502 / 0.1999 / 0.1507 / 0.1295 / 0.1000 / 0.0698 で、設定した 0.35 / 0.20 / 0.15 / 0.13 / 0.10 / 0.07 と一致しました。

結論。乱数生成器は rng.random() ひとつしかありません。4つの分布に共通の乱数源です。分布の情報は、乱数の側ではなく $F$ の形だけが持っています。 だから $F$ を書き下せる分布なら何でも作れる。「ライブラリに無い分布でも分析できる」という利点は、この一点に由来していました。

コラム：逆関数法と第4回の変数変換は同じ話なのか

第4回で変数変換とヤコビアンを扱いました。逆関数法もまた「一様乱数を別の値に変換する」操作です。同じ話なのか、ヤコビアンは出てくるのか。

同じ話です。逆関数法は変数変換の1次元・単調な場合の特殊ケースでした。そしてヤコビアンは出てきます。というより、ずっと使っていました。

第4回の公式は $f_X(x) = f_U(u) \cdot |du/dx|$ でした。逆関数法では $u = F(x)$ なので、

ヤコビアン $|du/dx|$ が、密度 $f(x)$ そのものです。 数値微分で確かめました。

さきほどの図で見た「累積分布関数の傾きが密度である」が、ヤコビアンでした。 1次元だと「傾き」と「密度」が同じものなので区別する必要がなく、ヤコビアンという言葉を持ち出さずに済んでいたわけです。2次元以上になると2つは別物になり、行列式を計算する必要が出てきます。それがこのあと出てくるボックス=ミュラー法の $r$ です。

棄却法：逆関数が書けないとき

指数分布は $F(x) = 1 - e^{-x}$ を $u$ について解けたので、1行で書けました。正規分布ではこれができません。

この積分は初等関数で書けず、誤差関数 $\mathrm{erf}$ が必要です。当然その逆関数 $F^{-1}$ も、$-\log(1-u)$ のような閉じた式にはなりません。

正直に書くと、正規分布は例外的に恵まれています。 数値ライブラリには erfinv があるので、実務では逆関数法も使えます。棄却法が本当に必要になるのは、自分で作った密度や正規化定数が分からない密度のときです。ベイズの事後分布がまさにそれで、第33回の MCMC（Markov Chain Monte Carlo＝マルコフ連鎖モンテカルロ法）につながります。ここでは仕組みが見やすいので正規分布を題材にします。

発想：密度を面積として見る

逆関数法は密度の情報を「累積分布関数の傾き」経由で使いました。棄却法はもっと素朴です。

1. 密度の曲線を、すっぽり覆う図形（まずは長方形）を用意する
2. その図形の中に一様にランダムな点をばらまく
3. 曲線より下に落ちた点だけ残し、その $x$ 座標を採用する

左は1400点をばらまいた例で、緑が採択435点、灰色が棄却965点です。採択された緑の点は、曲線が高いところに密集しています。 右は300万点でやって、採択した点の $x$ 座標だけのヒストグラムを描いたものです。

なぜ正しいのか。長方形の中に一様に点を打つと、どの領域に落ちる確率もその領域の面積に比例します。 だから $x$ の周りの細い縦の帯に落ちる確率は、その帯の中の「曲線の下の面積」に比例する。それはまさに $f(x) \times$ 帯の幅、つまり密度です。

採択率は「面積の比」そのもの

覆いが密度を完全に含んでいれば分子は 1 になるので、覆いの面積が小さいほど効率が良いことになります。

長方形 $[-4, 4]$ なら面積は $8 \times 0.398942 = 3.191538$ です。ただしこの場合、分子はぴったり1ではありません。$[-4,4]$ の外に $\pm 4$ より遠い裾がわずかに残っているからです。その分の質量を除くと 0.99993666 なので、

分子を 1 と置いて $1/3.191538 = 0.313329$ としてはいけません。 差は小数第5位ですが、この差の正体は「切り落とした裾」で、作られる分布が標準正規分布ではなくなっていることの現れです。次の表とジレンマの話に直結します。

長方形の幅を変えてみます。

ここにジレンマがあります。 広く取ると無駄が増えて採択率が落ちる。狭く取ると採択率は上がりますが、裾が切れて別の分布になってしまいます。 $[-1, 1]$ の採択率 0.856 は魅力的ですが、それで作られるのは「$\pm 1$ で打ち切られた正規分布」です。

正規分布の裾は無限に続くので、長方形では原理的に完全に覆えません。 どこかで切るしかない。

覆いを「密度に似た形」にする

長方形にこだわる必要はありません。覆いも確率分布でよいのです。ある定数 $c$ を使って

を満たす $g$ を用意します。すると採択率 $= 1/c$（$g$ が $f$ の全域を覆う確率分布である場合）。$c$ が 1 に近いほど良いということになります。

赤い余白がそのまま無駄になる分です。

長方形 $[-4,4]$ の行だけ $1/c$ と実際の採択率が食い違っています（0.313329 対 0.313309）。長方形は裾を覆いきれていないので $1/c$ が使えず、この差が切り落とした分です。$[-10,10]$ 以降は裾の質量が丸めの範囲で1になるため一致します。覆いが分布（コーシー・指数）なら全域を覆うので、採択率は厳密に $1/c$ です。

コーシー分布が使えるのは、裾が正規分布より厚いからです。 裾で追い越されないので全域を覆えます。「重い裾の分布で軽い裾の分布を覆う」は棄却法の定石で、逆はできません。

指数分布版（半正規分布を作って符号をランダムに付ける）で作った標準正規の検算は、平均 −0.00130、分散 0.99971、$P(|X| < 1.96) = 0.94998$（理論 0.95000）でした。コーシー版は平均 0.00027、分散 0.99920、$P(|X| < 1.96) = 0.95010$ です。

$c$ の求め方

$c$ は「$f/g$ の上限（sup）」です。厳密には最大値が達成されない場合もありますが、今回は $x=1$ で達成されます。
なお条件の本体は「$g$ の台が $f$ の台を含む」ことと「$\sup f/g < \infty$」であり、「重い裾で覆う」はそのための言い換えです（厳密には「$g$ の裾が $f$ より軽くないこと」で、同程度の減衰でも構いません）。そして $\int f = \int g = 1$ から必ず $c \ge 1$、つまり採択率は必ず1以下になります。

半正規分布を指数分布で覆う場合を実際にやってみます。

指数部 $-x^2/2 + x$ を微分して 0 にすると $-x + 1 = 0$、つまり $x = 1$ で最大です。

実測 0.760255 と一致しました。第1回から使っている「微分してゼロ」がここでも出てきます。

棄却法の限界：高次元で崩壊する

棄却法にも弱点があります。しかも、すでに何度も見た弱点です。

$d$ 次元の立方体 $[-1, 1]^d$ に一様に点を打ち、単位球の中に入ったものだけ採択する（＝球内の一様分布を作る）という問題を考えます。

$d = 15$ で200万点から19個しか採択されませんでした。

高次元では立方体のほとんどが「角」で、球はその中心のごく一部にしかありません。 覆いと密度の面積比が指数的に開くので、採択率が消えていきます。

これが第33回で MCMC が必要になった理由です。 高次元の分布から標本を得るのに、棄却法は使えません。MCMC は「棄却しても前の点に留まる」ことで、この崩壊を回避しています。逆に言えば、低次元では棄却法は今でも現役です。$d = 3$ なら1.91回に1回採択されるので、MCMC のような収束の心配がある道具を持ち出す理由がありません。適材適所でした。

ボックス=ミュラー法：無駄ゼロで正規乱数を2個作る

棄却法は無駄が出ます（長方形なら3.19回に1回）。逆関数法は書けません。そこで発想を変えます。1個ずつ作るのをやめて、2個セットで作る。

2次元標準正規分布（$x$ と $y$ が独立）の同時密度を書いてみます。

$x^2 + y^2$ しか出てきません。 これは「原点からの距離だけが効いていて、方向はまったく効いていない」という意味です。だから角度は完全に一様になります。

上の表の角度は arctan2 の $(-\pi, \pi]$ で測っているので平均が 0 になります。あとで出てくる手続きの $\theta = 2\pi u_2$ は $[0, 2\pi)$ で平均は $\pi$ ですが、幅が $2\pi$ の一様分布なら向きの取り方が違うだけで、標準偏差は同じ $\pi/\sqrt{3}$ です。

ここで話が終わっています。 分解して出てきた2つの部品は、どちらも簡単に作れるものでした。

• 角度 $\theta$ → 一様分布。$u$ をそのまま $2\pi$ 倍するだけ
• 距離の2乗 $S$ → 平均2の指数分布。逆関数法が使える（$S = -2\log u$）

正規分布そのものには逆関数法が使えないのに、2個セットにして極座標で見ると、逆関数法が使える部品に分解できたわけです。

なぜ極座標にすると指数分布が出るのか — ここがヤコビアン

第4回の変数変換では、変数を取り替えるときに「面積の伸び縮み」を補正する必要がありました。極座標では

この $r$ がヤコビアンです。ヤコビアン行列の行列式を数値で確認しておきます。

これを密度に掛けると、こうなります。

さらに $S = r^2$ と置くと $dS = 2r\,dr$ なので、

平均2の指数分布です。ヤコビアンの $r$ が、$S = r^2$ と置いたときの $2r\,dr$ にちょうど吸収されました。 これが「距離の2乗が指数分布になる」理由です。

なお $S = Z_1^2 + Z_2^2$ は定義から自由度2のカイ二乗分布です。つまり $\chi^2(2) =$ 平均2の指数分布。第7回で扱ったカイ二乗分布と、ここで出てきた指数分布が同じものだった、という同一視です。準1級では問われやすい形なので押さえておきます。

手続き

色つきの縦線が同心円に、灰色の横線が放射線に写ります。 $u_1$ が小さいほど大きな円になる（$-\log$ で反転するため）。「正方形を平面に巻きつける」変換だと見ると直感的です。

効率を比べます。

ヤコビアンの $r$ を落とすと壊れる

「$S = r^2$ を指数分布にする」を間違えて「$r$ 自体を指数分布にする」とどうなるか、実演しました。

注目してほしいのは分散の行です。誤り版でも分散 1.00147 で、ほぼ 1 になっています。 偶然ではなく、平均1の指数分布の2次モーメントが 2 なので $E[r^2] \cdot E[\cos^2\theta] = 2 \times 1/2 = 1$ になるためです。
この「分散だけ合う」現象は $r = -\log u$（平均1）にした場合に限ります。 $-2\log u$ をそのまま $r$ にすると $E[r^2] = 8$ で分散4になり、すぐに気づけます。
なお尖度も理論値が出ます。$E[r^4] = 24$、$E[\cos^4\theta] = 3/8$ なので $E[X^4] = 9$、つまり尖度は9。実測 9.1698 と整合しています。だから間違い方によって発見しやすさが変わるのが厄介なところです。

「平均と分散を確認したから合っている」では、この誤りを検出できません。 尖度（3 vs 9.17）と原点付近の割合（7.97% vs 21.83%）を見て初めて分かります。この連載で何度も踏んだ罠、理論値と一部が一致するとバグが隠れるが、ここにも出ていました。

擬似乱数は本当にランダムなのか

答えはまったくランダムではありません。 完全に決定的な漸化式です。線形合同法なら

同じ種（シード）を入れれば必ず同じ列が出ます。有限個の状態しかないので、いつか必ず一周します。

周期を実際に数えた

パラメータ次第で周期が $m$ よりずっと短くなります。 $a = 9, c = 0, m = 16$ は周期 2、つまり2つの値を往復するだけです。

結論として、現代の生成器では周期は問題になりません。 問題になるのは別の2つです。

本当の問題①：規則性は次元を上げると見える

RANDU の有名な欠陥を再現しました。この生成器は次の関係を厳密に満たしてしまいます。

実際に計算すると $9 - 6a + a^2 \bmod 2^{31} = 0$（$a = 65539$）でした。つまり3つ連続で取ると、必ず平面上に乗ります。

1次元の検定は全部通ります。 平均・分散・度数分布はどれも問題なし。欠陥は3次元にして初めて見えます。

これは連載で繰り返し出てきた構造と同じです。低次元の要約統計量が合っていても、高次元の構造は壊れていることがある。 そして今回のボックス=ミュラー法でも、$r$ を落とした版が「分散 1.00147」で合格してしまいました。同じ罠です。

本当の問題②：シードと再現性

落とし穴が2つあります。

「シードを固定したから再現できる」は、同じライブラリの同じバージョンでしか成り立ちません。 numpy も 1.17 で新しい API（default_rng ＝ PCG64）が入り、推奨が RandomState（メルセンヌツイスタ）から切り替わりました。後方互換のため np.random.rand などのトップレベル関数は今もメルセンヌツイスタなので、古いコードの出力が変わったわけではありません。数値を公開するときは生成器の名前も書くのが安全です。

もう1つは並行実行です。同じシードで複数プロセスを走らせると全部が同じ列を使い、実質的なサンプル数が増えません。numpy では SeedSequence や rng.spawn() で分岐させます。

第2部　乱数を「使う」技術

そもそも何のために乱数を使うのか

第1部で「任意の分布から乱数を作れる」ことが分かりました。しかしそれだけでは変な数字の列が手に入るだけです。そこから結論を出すには、もう1本の柱が必要でした。

この連載でずっとシミュレーションで検証してきた正当性は、②の大数の法則です。 ①だけでは入力が作れるだけ、②だけでは入力がない。両方が揃って初めて「シミュレーションで確かめた」と言えます。

そして①が必要になるのは、ライブラリに無い分布を扱うときです。実例で見ます。

題材：ブログ記事のPV

このブログの記事ごとのPVを題材にします。1記事あたりのPVはパレート分布（べき乗則）で表せることが多く知られています。最低10PV、形状パラメータ $\alpha = 1.5$ としました。

平均は中央値の1.89倍で、平均PVは半分以上の記事が到達しない値です。

見つかったこと1：標準偏差を報告してはいけない

「20本書いたら合計PVはいくらか」を平均 ± 標準偏差で出そうとしました。シードだけ変えて3回。

平均（596〜603）と中央値（489〜492）は安定しているのに、標準偏差だけが4倍以上ばらつきます。 $\alpha = 1.5$ のパレート分布は $\alpha < 2$ なので母分散が無限です。存在しない量を測ろうとしているので、値が落ち着きようがありません。

1回だけ回して「合計PVは約600、標準偏差は約2800」と書いていたら、次に測ったとき655になって説明できなくなります。シミュレーションを2回以上回すと、この種の事故が見つかります。

見つかったこと2：精度が $1/\sqrt{n}$ で上がらない

連載でずっと「精度はルートでしか上がらない」と扱ってきましたが、この分布では成り立ちません。標準偏差そのものが不安定なので、ばらつきは四分位範囲で測っています。

両対数での傾きは −0.318。$1/\sqrt{n}$ なら −0.5、$n^{-1/3}$ なら −0.3333 です。実測は後者に沿っています。

理屈はこうです。$\alpha < 2$ で母分散が無限なので、中心極限定理の前提（分散が有限）が崩れています。 $1 < \alpha < 2$ のとき、$\sum (X_i - \mu)$ の大きさは $\sqrt{n}$ ではなく $n^{1/\alpha}$ で伸びます（$\alpha$ 安定分布への収束。パレート分布では形状パラメータがそのまま裾指数になります）。だから標本平均のばらつきは $n^{1/\alpha - 1} = n^{-(1-1/\alpha)}$。$\alpha = 1.5$ なら $n^{-1/3}$ です。
$\alpha > 2$ なら分散が有限なので通常の $n^{-1/2}$ に戻り、$\alpha \le 1$ では平均そのものが存在しません（$\alpha = 2$ はちょうど $E[X^2]$ が発散する境界で、分散は無限のままです。sympy で確認しました）。 この式が意味を持つのは $1 < \alpha < 2$ のときだけです。

$n$ を10倍にしても、ばらつきは 3.16 分の1ではなく 2.15 分の1にしかなりません。 理論値は（$1 < \alpha < 2$ のもとで）$n^{-(1 - 1/\alpha)} = n^{-1/3}$ です。

この結論を理論だけで出すには安定分布と裾指数の理論が必要ですが、シミュレーションなら10行で出ます。そして実務上の意味は明確です。PVのような裾の重いデータでは、記事数を増やしても平均の精度が想定ほど上がりません。

見つかったこと3：ここで逆関数法が必須になる

上の分析には非現実的な前提があります。1記事が100万PVになる可能性を許していること。現実にはそんな記事は出ません。

そこで「上限 $b$ で切ったパレート分布」が必要になります。これは numpy に存在しない分布です。 逆関数法なら累積分布関数を $b$ で正規化して解き直すだけです。

上限を 1万 / 5万 / 100万 PV の3通りで比べます。（なおこのあとの節では、第1部の実演でも使った上限200PV版に戻します。 このブログの実測レンジに寄せた設定で、厳密値が手計算できるので検算に向いているためです。以降「上限200PV」と書いてある数値はすべてそちらです。）

中央値は 490.5 / 490.6 / 490.6 でほぼ完全に一致し、90%点も 825 / 827 / 827。ところが標準偏差は 412 → 684 → 1525 と3.7倍動きます。

上限をどこに置くかという検証しにくい仮定が、中央値には影響しないのに標準偏差にはまるごと乗ってしまう。だから中央値と90%点で語れば、仮定が結論に入りません。

これは「どの要約統計量を報告するか」という実務判断です。シミュレーションが答えを出したのではなく、仮定を1つ動かして、結論のどこが動くかを見せたわけです。

1変数なら数値積分に完敗する

ここで、私が最も思い違いをしていたことを書きます。モンテカルロ法は速くも正確でもありません。

「1記事の期待PV」を求めます。逆関数法で乱数を作って平均する方法と、同じ積分を等間隔グリッドで計算する方法を、同じ計算回数で比べます。

実は両者はまったく同じ積分を計算しています。

モンテカルロ法は $u$ を乱数で置き、グリッドは等間隔に置く。点の置き場所が違うだけです。 厳密値は 23.5551506580（上限200PVのパレート分布）。

100万回使って、グリッドは誤差 $5 \times 10^{-10}$、モンテカルロは $7 \times 10^{-2}$。1億倍以上の差でモンテカルロの負けです。

1つ補足します。上の表のモンテカルロ列は単一のシードでの実現値なので、たまたま大きく振れている可能性があります。この分布の標準偏差は厳密に $\sigma = 22.3288$（$E[X^2] = 1053.418$ から）なので、$n = 10^6$ での典型的な誤差は $\sigma/\sqrt{n} = 0.022329$。表の 0.0686 は3σ程度の外れです。典型値で比べても $0.022329 \div (4.91 \times 10^{-10}) \approx 4.5 \times 10^7$ 倍で、桁の主張は変わりません。

つまり1変数なら、乱数を使うのは点の置き場所をわざわざ雑にしているだけでした。等間隔に置いたほうがムラがなく、はるかに速く収束します。

逆転するのは次元

ではなぜモンテカルロ法が使われるのか。誤差の減り方の法則が違うからです。滑らかな被積分関数で、次元だけを変えて実測しました。

先に前提を書いておきます。ここでのグリッドは中点則（2次精度）のテンソル積で、$n^{-2/d}$ という指数はそのことと被積分関数が2階微分可能であることから出ます。指示関数のような不連続な量（確率の推定）では $n^{-1/d}$ に落ちるので、後で出てくる交差点の位置も変わります。

左（1次元）では青のグリッドが赤を圧倒し、右（8次元）では入れ替わっています。 要点は、赤の傾きが左右で変わっていないことです。

グリッドの指数は理論どおり $-2/d$ で、次元が上がるとどんどん 0 に近づきます（＝収束しなくなる）。モンテカルロは $-1/2$ のまま動きません。

同じ計算回数（おおむね $n \approx 10^6$）での誤差を並べます。グリッドは軸あたりの点数 $p$ に対して $n = p^d$ しか取れないので、$n$ をぴったり揃えられません。 $d=1$ は $n = 1869$（1次元では点を増やすと誤差が浮動小数の下限に達するので打ち切り）、$d=2$ は 1,144,900、$d=4$ は 1,185,921、$d=8$ は 390,625（$= 5^8$）です。

この滑らかな例で交差するのは $d = 4$ です。 不連続な量（確率）ならグリッドは $n^{-1/d}$ なので、交差点は $d = 2$ まで下がります。

次元が上がるとグリッドは軸ごとに点を掛け算しなければならず、点が指数関数的に必要になります。$d = 20$ で軸10点なら $10^{20}$ 回の評価で、1秒に10億回計算しても3175年かかります。同じ問いにモンテカルロは200万回・数秒で答えます。

これが「$1/\sqrt{n}$ が高次元で有利になる」の中身でした。$1/\sqrt{n}$ は「遅い」のではなく「次元に鈍感」という長所だったのです。

モンテカルロ法も次元の呪いを受ける（ただし別の場所で）

では、モンテカルロ法は次元の呪いから自由なのか。自由ではありません。ただし呪いが出る場所が違います。

モンテカルロ法の誤差はこう分解できます。

グリッドと比べたのは $1/\sqrt{n}$ の部分で、これは次元に対して不変でした。呪いは $\sigma$ のほうに出ます。 そして $\sigma$ が次元とともにどう動くかは、何を推定しようとしているかで正反対になります。

3つとも閉じた式が出るので、厳密値を並べます。A は $1/\sqrt{d}$、B は $\sqrt{(\pi^2/8)^d - 1}$、C は $\sqrt{2^d - 1}$ です。

A の実測値（$n = 5$ 万）は 0.9987 / 0.7109 / 0.5031 / 0.3553 / 0.2506 / 0.1769 で厳密値と一致しました。ところが B を標本標準偏差で推定すると、$d = 32$ だけ 42.83 という値が出ます。厳密値 28.78 の1.5倍です。

シードを変えて8回測り直したら 42.83 / 22.41 / 21.19 / 28.64 / 76.24 / 26.34 / 22.33 / 33.18 と3倍以上ばらけました（$d = 16$ では 5.19〜5.33 で安定していて、厳密値 5.272 の近傍にいます）。原因は $\prod \cos$ の $d = 32$ での尖度が $(3/2)^{32} = 4.3 \times 10^5$ もあることで、$\sigma$ そのものが推定できない領域に入っています。

この表を作るとき、私自身がこの章の主題を踏みました。 「$\sigma$ が指数的に増える」という定性的な主張は正しいのですが、その $\sigma$ を標本標準偏差で1回測って表に載せていました。母分散が無限のときに標準偏差が暴れる話（さきほどのパレート分布）と、まったく同じ構造です。閉じた式があるなら使うべきでした。

A：次元が上がるほど良くなるケース

「$d$ 個の値の平均」を推定する問いです（ここでは1個を平均1の指数分布としたモデル例で測っています）。1回の試行で $d$ 個の平均を取るので、試行の中身自体が $d$ 個の平均化になっており（独立なら）$\sigma$ が $1/\sqrt{d}$ で縮みます。 実測でも $d=32$ で $0.9987 \div \sqrt{32} = 0.1765$ に対して 0.1769 でした。

ただし正確に言うと、A は「同じ積分の次元だけを変えた」のではなく、$d$ が上がると推定したい量そのものが安定する型です。「高次元の積分が楽になる」のではありません。

ここまで見てきた20本の記事の問い（合計・分位点・シェア）は、すべてこのA型でした。 だから20次元でも40万回で安定した答えが出ていたわけです。

B：積の形はまずい

$d$ 個の量の積の期待値 $E\left[\prod_i h(X_i)\right]$ を推定する型です。ここでは $\int_{[0,1]^d} \prod_i \cos(\pi u_i / 2)\, du$（真値 $(2/\pi)^d$）を使いました。

積は1つでも小さい値を引くと全体が小さくなるので、1試行の値が桁で暴れます。 そのうえ答え $I$ 自体が $d$ とともに指数的に小さくなるため、$\sigma \div I = \sqrt{(\pi^2/8)^d - 1}$ が指数的に増えていきます（$\pi^2/8 = 1.2337$ なので、$d$ が1増えるごとに約1.11倍）。

実務でこの形が出るのは尤度の積や経路の確率です。対策は「対数を取って和に直す」か「逐次モンテカルロ法（粒子フィルタ）で途中でリサンプリングする」です。

C：次元が上がるほど悪くなるケース

「$d$ 次元の一様乱数がすべて 0.5 未満」という確率 $0.5^d$ を素朴に数えます。

$d = 20$ で命中が1件も出ず、推定値が 0 になります。

$d = 16$ の相対誤差 1.70% が $d = 12$ の 5.68% より良く見えますが、これは偶然です。期待される命中数が 15.26 件のところ実際に 15 件出ただけで、次に回せば 8 件や 22 件になります。命中数が2桁を切ったら、相対誤差の数値は運です。

相対誤差10%に必要な $n \approx 100/p$ で計算すると、

指数関数的な悪化、つまりまさに次元の呪いです。

それでも「傾き」は3つとも $-1/2$

ここが精密にしておきたい点です。A も B も C も、$n$ を増やしたときの減り方の法則は同じでした。250回反復して二乗平均平方根誤差（RMSE = Root Mean Squared Error）を測りました。

3つとも傾きは $-0.50$ で、灰色の $1/\sqrt{n}$ の線と平行です。違うのは線の高さ（出発点）だけでした。

言い換えるとこうなります。

モンテカルロ法の呪いは「$n$ を増やせば必ず解決する。ただし必要な $n$ が天文学的」という形をとります。 体質は健全なのに要求量が非現実的。グリッドは $n$ を増やしても体質自体が治りません。

20次元の問いに、何が答えられるか

「$F$ が計算できるなら、わざわざ乱数を使う必要はないのでは」という疑問が自然に出ます。半分は正しく、上で見たとおり1変数なら乱数は不要です。

分かれ目は「$F$ が計算できるか」ではなく「知りたい量が何変数の積分か」でした。

同じ設定（上限200PVのパレート分布・20本）で、式で書ける問いと書けない問いを並べます。1記事の期待PVは $E[X] = 23.5552$（厳密）です。

式で書ける問い（乱数は不要）

期待値は線形なので次元が実質1に落ちます。最大値も「全部が150以下」に分解できて積になります。

$F$ が完全に分かっていても式が書けない問い

この2つの表の対比が答えです。 合計PVの期待値は1行の式で出ます（471.103）。しかし合計PVの中央値（456.4）は出ません。期待値は足し算に分解できるが、分位点は分布の形そのものを要求するからです。

そして実務で使うのは分位点のほうでした。期待値 471.103 は中央値 456.4 より大きく、半分以上の「20本セット」が到達しない値だからです。

本数を変えたときの見通し（上限200PVのパレート分布・20万回。すべて20〜80次元の積分で、1つも式では出ません。20本の行が上の表と 456.4 対 456.2 でわずかに違うのは、反復回数とシードが違うためです）

本数を8倍にすると合計の中央値は 8.5 倍（ほぼ比例）ですが、変動係数は 0.301 → 0.106 と3分の1に落ち、1本の当たり記事への依存度も 24.7% → 7.2% に下がります。 「本数を増やすと当たり待ちのギャンブル性が薄まる」という定量的な結論です。

$F$ が計算できることは、逆関数法を使わない理由ではなく逆関数法を使える前提のほうでした。$F$ が書けないと逆関数法自体が動きません。

分散減少法：乱数の置き方を変える

$1/\sqrt{n}$ は、母分散が有限で、点を独立に $n$ 個引くかぎり変えられません。この2つの条件が破れると話が変わります。母分散が無限なら遅くなる（さきほどのパレート分布の $n^{-1/3}$）。そして点を独立に引くのをやめれば速くなります。

つまり分散減少法には2つの型があります。

それが分散減少法という一群の技術でした。

素朴なモンテカルロ法の弱点は、乱数がムラを作ることです。24個の点を置いた図を並べます。

素朴（赤）は目で見て分かるほどムラがあります。 0.05〜0.15 が空白なのに 0.2 付近と 0.6 付近に固まっている。下に行くほどムラが消え、最下段のグリッドは完全に均等です。分散減少法はこの赤を青に近づける工夫だと考えると全体が見通せます。

期待値を推定するかぎり、4つとも推定値に偏りを作りません（重点サンプリングも重みで正確に補正します）。ばらつきだけを削るので、「タダで精度が上がる」ように見えます。

ただし条件が3つ付きます。重点サンプリングは提案分布の台が元の分布の台を覆うこと（後述の「6だけ出るサイコロ」は目標が $\{6\}$ だから救われているだけです）。制御変量法は係数 $b$ を別の標本で決めること（同じ標本から推定すると小さな偏りが出ます）。そして中央値や分位点の推定は、素朴法を含めてどの方法でも有限標本では厳密に無偏ではありません。

実測：効くが、問いによって桁が違う

同じ問い（20本の記事の合計PV、上限200PVのパレート分布）で、目標だけを変えて比べました。$n = 100{,}000$ を120回反復し、推定値の標準偏差を測っています。

目標が「和の期待値」のとき（線形）——厳密値は 471.1030 です。

層別サンプリングが9707倍。 ただしこれは $n = 100{,}000$ という一点での倍率です。ここでの層別サンプリングは、20次元それぞれを独立に $n$ 分割してシャッフルするラテン超方格法として実装しています。$n$ を変えて測ると倍率自体が動きます。

（この表は40反復での測定です。上の表は120反復なので $n=10^5$ の値がわずかに違います。）

倍率が $n$ にほぼ比例して開いていきます。 両対数の傾きを測ると素朴が −0.531（＝$n^{-1/2}$）に対して層別は −1.477（＝$n^{-3/2}$）でした。これは定数倍の改善ではなく収束の次数そのものが変わっているということです。

なぜこの目標でだけ桁違いに効くのか。推定している量 $S = \sum_k F^{-1}(u_k)$ が完全に加法的（各座標の寄与を足しただけ）だからです。ラテン超方格法は座標ごとのムラを消すので、加法的な成分の分散をほぼ丸ごと落とせます。

目標が「和の中央値」のとき（非線形）

目標が「90%点」「600超の確率」のとき

同じ手法・同じデータなのに、改善倍率が 9707 倍から 1.75 倍へ落ちました。

層別サンプリングが劇的に効いたのは、目標が期待値（＝足し算）だったからです。各区画の寄与を足すという構造が層別と完全に噛み合う。ところが中央値や分位点は「並べ替えて真ん中を取る」操作で、区画ごとに分解できません。分散減少法の効果は、推定したい量の構造に依存します。

そして実務で使うのは中央値や90%点のほうでした。つまり「1.75倍のために実装を複雑にするか」という判断になります。

重点サンプリング：わざと歪めて重みで割り戻す

分散減少法の中で、重点サンプリングだけは別格です。上の4つの中では、まれな事象に効く唯一の手段だからです（ほかに分割法や交差エントロピー法がありますが、いずれも「事象の近くを重点的に引く」という同じ発想です）。ただし理屈が入りにくいので、サイコロ1個から始めます。

重みとは何か

歪んだサイコロがあって、6が出る確率を推定したいとします。素朴にやると100回振って平均2回しか6が出ないので精度が出ません。そこで6が出やすいサイコロに作り替えて振ります。

重みの意味はこれだけです。

本来 0.02 しか出ないはずの目を 0.50 で出させたので、25倍ズルをした。だから1回を $1/25 = 0.04$ 回分として数える。

選挙の出口調査で、ある投票所だけ意図的に多めに調査したら、集計のときに人数比で割り戻すのと同じです。第21回の「抽出確率が違うサンプルは重みをつけて集計する」がそのまま出てきています。

提案のサイコロを10回振って、「6が出たときだけ 0.04 を記録、それ以外は 0」を記録します。

平均 $= 0.04 \times 5 \div 10 = 0.0200$。真値 0.02 と一致しました（この10回はちょうど6が5回出た当たり回です）。

なぜ合うのか。提案では6が10回中およそ5回（$= 0.50 \times 10$）出ます。その5回に 0.04 を掛けると $0.50 \times 0.04 = 0.02$。$q_6$ と $w_6 = p_6/q_6$ を掛けると $q_6$ が約分されて $p_6$ だけが残る。 これが「重みで補正する」の中身です。

どちらも 0.02 に収束し（偏りなし）、標準誤差は7倍改善しています。

理想の提案分布は作れない

「6を出やすく」がうまくいったなら、もっと極端にしたらどうなるか。6しか出ないサイコロまで行ってみます。

6しか出ないサイコロを使うと、毎回必ず 0.02 が記録されます。全部同じ値なのでばらつきがゼロ。1回振るだけで厳密に答えが出ます。

ところがこれは使えません。重み $w_6 = 0.02$ を計算するには $p_6 = 0.02$ を知っている必要があり、そして $p_6$ こそが求めたい答えでした。

理想の提案分布は「答えを知っていれば、答えがそのまま重みとして毎回出てくる」というものです。すでに知っている数字を読み上げているだけなので、実際には作れません。

事象が複数の目にまたがる場合で書くと構造が見えます。「5か6が出る確率」（真値 0.10）なら、理想の提案分布は次の3ステップで作れます。

1. 元の分布から事象の外（1〜4）を捨てる → 残るのは 0.08 と 0.02。合計 0.10 しかないので確率分布になっていない
2. 0.10 で割って合計1に直す → 0.8 と 0.2。これが「$\{5,6\}$ に条件づけた分布」
3. この分布を提案に使うと、重みは $0.08 \div 0.8 = 0.10$、$0.02 \div 0.2 = 0.10$。どちらも 0.10 で一定

②で割った数がそのまま重みとして戻ってくるだけなので、どの目でも同じ値になります。教科書の $w = f(x) / (f(x)/p) = p$ はこれを1行で書いたもので、$p$ は「割るときに使った 0.10」つまり求めたい答えそのものです。分散は素朴法の 9.0000e−02 から 0 になります。

重みは答えを知らなくても計算できる

ここで混乱しやすい点を1つ潰しておきます。「$p$ が分からないなら重みも決められないのでは」という疑問です。

上のサイコロの例が悪かったので補足します。あの例では「求めたい答え $P(X=6)$」と「重みの計算に必要な個々の確率 $p_6$」が同じ数字でした。だから循環しているように見える。しかも答えが仕様表に書いてあるので、そもそもシミュレーションする必要がありませんでした。

サイコロを3個に増やすと分離します。

問いは「このサイコロを3個振って、合計が16以上になる確率」。216通りを全列挙した厳密値は 0.00060800 で、これは仕様表のどこにも書いてありません。提案は大きい目が出やすいサイコロ $q = (0.05, 0.05, 0.10, 0.15, 0.25, 0.40)$ にします。

この表の全部の数字は、仕様表の12個（$p$ の6個 + $q$ の6個）から作られています。0.000608 という答えは1度も使っていません。

素朴法は $n = 10{,}000$ で命中3回しかなく、推定値が真値の半分になっています。

一般的な構図はこうです。重点サンプリングが要求するのは「観測した1点に元の分布が割り当てる確率（密度）を計算できること」だけ。要求しないのは「その分布から導かれる複雑な量（裾確率・分位点・和の分布）を知っていること」。 式に代入できる量と、式から導かねばならない量の差でした。

そして循環が起きる設定は、シミュレーションが不要な設定と同じでした。

重みが元の世界を厳密に復元することの確認

コイン10回で「表が9回以上出る確率」（真値 $11/1024 = 0.0107421875$）を題材にすると、重みの働きが厳密に確認できます。表の確率 0.8 のイカサマコインを提案に使い、表が $k$ 回出た列に対して

各 $k$ について「イカサマでの確率 × 重み」を計算すると、公平コインでの確率にぴったり戻ります。

11個すべてで一致しました（差はすべて浮動小数の丸め誤差）。重みを掛ける操作は、イカサマコインの世界を公平コインの世界に厳密に戻す変換になっています。 だから足し合わせれば答えが出る。$n = 10^7$ での推定は重点 0.01074752（標準誤差 4.99e−06・命中率 37.6%）、素朴 0.01070720（標準誤差 3.25e−05・命中率 1.1%）でした。

まれな事象での威力

いよいよ本題です。「Exp(1) を20個足した合計が 100 を超える確率」を推定します。ガンマ分布の上側確率なので厳密値が出ます：3.764894e−23。

素朴法は1000万回で命中0件。重点サンプリングは1000回で桁が合い、10万回で誤差1%。$10^{-23}$ という確率を10万回のシミュレーションで測っています。 命中率は 0% から 47% に上がりました。

目安 $\theta =$ しきい値 $\div\, d$ はどこから来るのか

上の $\theta = 5$ は当てずっぽうではありません。$100 \div 20 = 5$ です。この目安の出どころを確認します。しきい値 $a$、変数の個数 $d$ とすると $\theta^* = a/d$ になります。

まず理想の分布の姿を測る

理想の提案分布（事象に条件づけた分布）は作れませんが、どんな形をしているかは測れます。$S > 40$ の場合で、棄却法を使って条件を満たす標本だけを4億回引いて集めました。

「合計が40を超えた世界」では、1本あたりの平均が 1 ではなく 2 前後になっています。 しきい値をちょうど超えるあたりに集まるので $a/d$。つまり $\theta = a/d$ は「理想の分布の平均を真似た設定」でした。形まで一致させることはできませんが、中心は合わせられます。

母関数が出てくる

なぜ「平均 $\theta$ の指数分布」で真似られるのか。ここで第9回の母関数が効きます。元の密度に $e^{tx}$ を掛けて正規化する操作を指数傾斜といいます。

これは平均 $1/(1-t)$ の指数分布です。つまり「平均 $\theta$ の指数分布から引く」という操作は、実は $t = 1 - 1/\theta$ で傾斜させたことと同じでした。指数分布は傾斜させても指数分布のまま（指数型分布族の性質）なので、族の中で閉じています。

そして重みを書き直すと、モーメント母関数 $M(t) = 1/(1-t)$ がそのまま現れます。

密度比を素直に計算した値とこの式の値を比べたところ、最大の差 8.327e−16（浮動小数の丸め誤差）で一致しました。

最適な $t$ は、傾斜した世界での $S$ の期待値をしきい値に等しくする条件で決まります。

日本語に直すと「まれな事象を、ふつうの出来事に変える」です。元の世界では合計40は平均20から遠い外れ値ですが、$\theta = 2$ の世界では合計40が平均そのものになります。

実測した $S$ の平均は、元の分布で 20.01、$\theta = 2$ で 40.00、$\theta = 6$ で 119.99 でした。$\theta = 6$ まで行き過ぎると、しきい値40の近傍が逆に空きます。 行き過ぎが悪いのは「まれな事象を今度は反対側からまれにしてしまう」ためです。