JA EN
体系確率・統計
·★ 会員·12分で読めます

モンテカルロ法を1から — 積分をサイコロで解く

解けない積分は、乱数を降らせて平均を取れば近似できます。大数の法則がなぜ 1/√N の遅さを生むのか、重点サンプリングで何が救えるのか、そして正規化定数が分からない分布をMCMCがどう扱うのかを、前提知識ゼロから積み上げます。

対象textタスクmath

池の面積は、雨粒を数えれば分かる

100平方メートルの敷地の真ん中に、いびつな形の池があります。この池の面積を、測量せずに知りたい。

方法があります。敷地全体に均等に雨が降るのを待って、落ちた雨粒のうち何割が池に落ちたかを数えるのです。1000粒のうち283粒が池なら、面積はおよそ28.3平方メートル。池の形を数式で書き下す必要も、複雑な図形の公式を思い出す必要もありません。「当たった割合」を「大きさ」に読み替える、それだけです。

これがモンテカルロ法の全部です。1940年代、ロスアラモスの数学者スタニスワフ・ウラムが病床でソリティアの勝率を考えていて、「場合の数を数え上げるより、実際に何百回も配って数えたほうが早い」と気づいたのが発端でした。フォン・ノイマンとともに中性子の拡散計算に応用し、モナコのカジノ地区にちなんで「モンテカルロ」と名付けられています。賭場の名前が付いた計算法ですが、やっていることは雨粒を数えるのと変わりません。

積分は「平均」だと読み替える

学校では、積分は「細長い短冊を敷き詰めて足す」と習います。モンテカルロ法はここを別の角度から見ます。積分は、幅 × 平均の高さである、と。

abf(x)dx=(ba)E[f(U)],UUniform(a,b)\int_a^b f(x)\,dx = (b-a)\cdot \mathbb{E}\big[f(U)\big], \qquad U \sim \mathrm{Uniform}(a,b)
(1)

記号を1つずつ読みます。E\mathbb{E} は「平均」、UUniform(a,b)U \sim \mathrm{Uniform}(a,b) は「aa から bb の間で、どこも同じ確からしさで選んだ1点」、f(U)f(U) はそこでの関数の高さです。つまりこの式は、区間の中からでたらめに1点選んだときの高さの平均に、区間の幅を掛ければ面積になると言っているだけです。長方形の面積が「幅 × 高さ」なのと同じ理屈で、高さが場所によって変わるぶんを平均でならしています。

平均は、実際に取ってきた点で代用できます。

I^N=baNi=1Nf(xi)\hat{I}_N = \frac{b-a}{N}\sum_{i=1}^{N} f(x_i)
(2)

x1,,xNx_1,\dots,x_N を区間からでたらめに NN 個選び、そこでの高さを足して個数で割る。要するに「何点かで味見して、その平均を全体の代表とみなす」という、選挙の出口調査と同じ発想です。これがモンテカルロ積分で、ff がどんなに複雑でも、値さえ計算できれば使えます。

なぜ短冊を敷き詰めないのか

1次元なら、短冊を等間隔に並べる(台形則やシンプソン則)ほうが正確です。わざわざ乱数を使う理由はありません。事情が変わるのは次元が増えたときです。

1軸あたり10点の格子を作るとします。2次元なら 102=10010^2 = 100 点、3次元なら1000点。ここまでは平気です。しかし10次元なら100億点、20次元なら 102010^{20} 点。1軸を2分割する最低限の粗さでも、2202^{20} で100万点を超えます。格子の点数は次元数に対して指数関数的に増える — これが次元の呪いです。

FIG 1横軸の n を「次元の数」と読み替えてください。1軸を2分割するだけの最も粗い格子でも点の数は O(2ⁿ) で、20次元あたりで他の曲線は全部床に貼り付きます。縦軸を線形に切り替えると、この差が「少し多い」ではなく「桁が違う」ことが見えます

モンテカルロ法の強みはここにあります。あとで見るように、この方法の誤差は 1/N1/\sqrt{N} で減り、その速さは次元数に一切依存しません。100次元だろうと1万次元だろうと、1万点なら1万点なりの精度が出ます。1次元では格子に負ける手法が、高次元では唯一現実的な手段になる。金融のポートフォリオ評価も、CGのレンダリングも、ベイズ推定も、扱う積分が軒並み高次元だからモンテカルロ法を使っているのです。

1/√N の正体 — 大数の法則と中心極限定理

「たくさん取れば真の値に近づく」と保証するのが大数の法則、「どれくらいの速さで近づくのか」を答えるのが中心極限定理です。実務で効くのは後者です。

I^N\hat{I}_N は独立な NN 個の値の平均なので、その分散は元の分散の NN 分の1になります。

Var(I^N)=σ2N標準誤差=σN\mathrm{Var}(\hat{I}_N) = \frac{\sigma^2}{N} \quad\Longrightarrow\quad \text{標準誤差} = \frac{\sigma}{\sqrt{N}}
(3)

σ\sigmaf(U)f(U) のばらつき(標準偏差)です。分散が NN で割られても、誤差として感じるのはその平方根なので、誤差は N\sqrt{N} でしか小さくならない。ここが肝心です。

数字にすると痛みが分かります。円周率をダーツで求める例(正方形に点を打ち、四分円に入った割合を4倍する)では、標準誤差がおよそ 。小数第2位を安定させるには 万回、もう1桁欲しければ約270万回です。精度を1桁上げるコストは10倍ではなく100倍。

この先にあるもの

§

ここから先は会員限定です

解説記事371本・教科書26章・学生モード48単元・論文精読6本が、月額¥490ですべて読み放題になります。新しい解説は毎日3本ずつ増えます。いつでも解約でき、解約後も期間の終わりまで読めます。

会員の方はログインすると続きが表示されます

コメント

コメントにはログインが必要です