モンテカルロ法のイメージをアニメーションで掴む

2020/07/20

昔モンテカルロ法を勉強していた頃のノートを発掘したのでメモします。 ここでは厳密な議論は避けてモンテカルロ法のざっくりとしたイメージの紹介に専念します。

モンテカルロ法とは

Wikipedia を見ると、

シミュレーションや数値計算を乱数を用いて行う手法の総称

とあります。有名なのはダーツを投げて円周率を求める計算ですよね。とりあえずやってみましょう。

1×1 の領域で一様に分布する乱数をプロットしていき、1/4 円の中に入ったプロットを数え上げていきます。この時プロット数が十分に大きければ

(円の面積)(正方形の面積)(円内のプロット数)(全プロット数)\small \frac{(\text{円の面積})}{(\text{正方形の面積})} \sim \frac{(\text{円内のプロット数})}{(\text{全プロット数})}

になることを利用して円周率(要は 1/4 円の面積)の計算をしています。

面積を計算する

上記では乱数を用いて面積比を求めていましたが、今度は 1/4 円の面積をより直接的に計算してみます。

区間[0,1]で一様に分布する乱数X1,X2,...,XNX_1,X_2,...,X_Nを考えてみます。この時サンプル数NNが十分に大きければ1/4 円の面積SSは、1X2\sqrt{1-X^2}の値を足し合わせた値で書けそうです。すなわち、

S=011x2dx1N{1X12+...+1XN2)}\begin{aligned} S &= \int^1_0 \sqrt{1-x^2} dx \\ &\sim \frac{1}{N} {\small \left\{\sqrt{1-X_1^2}+...+\sqrt{1-X_N^2}) \right\}} \end{aligned}

で書けそうな感じがします。実際に計算してみると、

上記ではπ(=4S)\pi (=4S)を表示していますが、うまく計算できていそうなことが分かると思います。

一般的に関数f(x)f(x)の区間[a,b]での積分は、区間 [a, b]の一様分布 pr(x)p_r(x)を導入すると、

abf(x)dx=abf(x)pr(x)pr(x)dx1Ni=1Nf(Xi)pr(Xi)=(ba)Ni=1Nf(Xi)\begin{aligned} \int^b_a f(x) dx &= \int^b_a \frac{f(x)}{p_r(x)} p_r(x)dx \\ &\sim \frac{1}{N} \sum^N_{i=1} \frac{f(X_i)}{p_r(X_i)} \\ &= \frac{(b-a)}{N} \sum^N_{i=1} f(X_i) \end{aligned}

で書くことができます。ここで、X1,X2,...,XNX_1,X_2,...,X_Nは一様分布pr(x)p_r(x)からサンプリング(確率分布からサンプルを得ること)してきた値で、また区間 [a, b]の一様分布pr(x)=1/(ba)p_r(x)=1/(b-a)であることを利用しています。

このような乱数を用いた数値積分はモンテカルロ積分と呼ばれ数値計算の分野ではよく利用されます。

ところで

上記の円周率の計算はたいして精度よくないです。これはサンプル数NNが十分ではないためです。モンテカルロ積分の誤差はサンプル数NNに対してO(1/N)O(1/\sqrt{N})で書けます。実際のところ、上記のような低次元の積分ではモンテカルロ積分よりも台形公式に代表される決定論的な手法の方が効率が良いです。

一方で、高次元(Davis と Rabinowitz によれば 15 次元以上)の積分を行う場合、モンテカルロ積分は強力な手法になります。これは原理的にモンテカルロ積分の誤差が非積分関数の次元に依存しないためです。台形公式などの求積法ではいわゆる次元の呪いの影響を強く受けてしまいます。

期待値を計算する、あるいはサンプリング

このモンテカルロ積分はよく期待値IIの計算に利用されます。すなわち、

I=f(x)p(x)dx1Ni=1Nf(Xi)I = \int f(x)p(x) dx \sim \frac{1}{N} \sum^N_{i=1} f(X_i)

です。つまりモンテカルロ積分の考えによれば、任意の確率分布p(x)p(x)に従う十分な数のサンプルX1,X2,...,XNX_1,X_2,...,X_Nを得ることができればは平均するだけでその期待値が得られるわけです。

ではどうすれば任意の確率分布からサンプリングが行えるでしょうか。

一様乱数

今更ですが一様乱数にも触れておきます。一様乱数でヒストグラムを作ってみると、

となり一様分布からサンプリングできていることが分かります(もちろんサンプル数を増やせばより一様になる)。

一様乱数の発生方法は、線形合同法から Mersenne TwisterXorshift など様々ありますがここでは詳細には触れません(結構沼なので...)。今回は文明の利器としてありがたく利用させて頂くことにしましょう。

逆関数法

逆関数法は、累積分布関数の逆関数を用いるサンプリング手法です。目的とする確率分布p(x)p(x)の累積分布F(x)F(x)を用いるとサンプルXXは、

X=F1(U)X = F^{-1}(U)

でかけることを利用します。ここでUUは区間[0,1)の一様乱数です。

例えば、指数分布の累積分布の逆関数は、

F1(x)=log(1x)λF^{-1}(x) = -\frac{\log{(1-x)}}{\lambda}

で書けます。これを利用してヒストグラムを作ってみると、

となり指数分布からサンプリングできていることが分かります。

棄却サンプリング

逆関数法ではその名の通り逆関数を利用しましたが、解析的に逆関数が得られないような少々複雑な分布(正規化定数が分からない分布とかも)からサンプリングしたい場合も多々あります。棄却サンプリングでは、このようなサンプリングを提案分布を導入することで実現します。

目的とする確率分布p(x)=p~(x)Zp(x)= \frac{\tilde{p}(x)}{Z}を覆うような提案分布q(x)q(x)を導入します。すなわち、

p~(x)<kq(x)\tilde{p}(x) < kq(x)

を満たすq(x)q(x)です。ここで、k は任意の定数です。このとき棄却サンプリングは、

  1. 提案分布q(x)q(x)に従う乱数XXを生成する。
  2. 一様乱数YY0Y<q(X)0 \leq Y < q(X)で発生させる
  3. Y<p~(X)Y<\tilde{p}(X)であればXXを確率分布p(x)p(x)のサンプルとして採択し、そうでなければ棄却する。

という手順で進みます。

例えば、分布p~(x)=(1x2)\tilde{p}(x)=\sqrt{(1-x^2)}を一様分布から棄却サンプリングするとことを考えると、

  1. 一様乱数u1,u2u_1, u_2を発生させる。
  2. u1<p~(u2)u_1<\tilde{p}(u_2)であればu2u_2を確率分布p(x)p(x)のサンプルとして採択し、そうでなければ棄却する。

となります。わお!これはあの面積比から円周率を求める計算と同じですね!そうですあの円周率の計算は棄却サンプリングの一種なんです。ここから得られるヒストグラムを作ってみると、

となります。だいたい目的とする分布からサンプリングできていることが分かります。

マルコフ連鎖モンテカルロ法

棄却サンプリングでは提案分布を導入することで、任意の確率分布からのサンプリングを実現しました。しかし、高次元な分布など非常に複雑な分布に対しては適切な提案分布を選ぶのは困難(というかほぼ不可能)になります。そんなときに使うのがMCMC(マルコフ連鎖モンテカルロ法)です。

MCMC は任意の確率分布p(x)p(x)定常分布に持つマルコフ連鎖を生成することで、確率分布p(x)p(x)からサンプリングする手法です。かいつまんで言ってしまうと、目的とする確率分布に従って状態が生成されるように遷移確率を定めて、この遷移確率のもと状態を次々と遷移させていくことでサンプルを得る手法です。

詳細つり合い条件

遷移確率の決め方として一般的に利用されるのが詳細つり合い条件です。詳細つり合い条件は、ある確率分布p(x)p(x)における状態AAから状態BBへの遷移確率をP(AB)P(A\rightarrow B)とおくと、

pAP(AB)=pBP(BA)p_A P(A\rightarrow B) = p_B P(B\rightarrow A)

でかけます。この条件のイメージはこの記事の説明が分かりやすくて、

pA,pBp_A, p_Bはそれぞれの国にいる人の人口で、P(AB),P(BA)P(A\rightarrow B), P(B\rightarrow A)は、A 国にいる人のうちどのくらいの割合が B 国に移住するかを表している。上記の条件が満たされていれば、移住(遷移)の後でも人口が変わらないのがわかると思う。この条件を詳細釣り合い条件と呼ぶ。

ここで重要なことは遷移によって各状態の存在確率pA,pBp_A, p_Bが変化していないことです。 この詳細つり合い条件を満たすように遷移をくり繰り返すことで目的の定常分布が得られることが知られています(詳細な証明は割愛)。

メトロポリス法

詳細つり合い条件から遷移確率を得る方法はいくつかありますが、ここではシンプルなメトロポリス法を紹介します。

メトロポリス法は、

  1. 状態xxを与える。
  2. 一様乱数等を用いて新たな状態x+dxx+dxを生成する。
  3. 遷移確率の評価をし採択ならx+dxx+dxを新たな状態とし、棄却ならxxを新たな状態とする。

という手順で進みます。この時、遷移確率は以下を用います。

P(AB)=min{1,pBpA} P(A\rightarrow B) = \min{ \left\{ 1, \frac{p_B}{p_A} \right\}}

例えば、標準正規分布の遷移確率は、

P(xx)=min{1,e(x2x2)/2}P(x\rightarrow x') = \min{ \left\{ 1, e^{(x^2-x'^2)/2} \right\}}

になります。これから得られるヒストグラムを作ってみると、

となります。大方のところ標準正規分布からサンプリングできていることが分かります。ただし、MCMC ではサンプル間に相関があることには注意が必要です。

おわりに

使用したプログラムは(雑ですが)下記にあげています。


必要であれば参考にしてみてください。(.....誤りなどがもしあれば PR 等いただけるととても嬉しいです。)

参考

logo

たまおさ

釣りとか登山とか好きです。(@tamaki_osamu)