モンテカルロ法

公開

目次

モンテカルロ法とは

Wikipediaより、

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

とあります。有名なのはダーツを投げて円周率を求める計算。

mc-dart

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

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

になることを利用して円周率(要は 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}

で書けそう。実際に計算してみると、

mc-int

上記ではπ(=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を得ることができればは平均するだけでその期待値が得られる

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

一様乱数

一様乱数でヒストグラムを作ってみると、

mc-uniform

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

一様乱数の発生方法は、線形合同法から 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}

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

mc-inverse

棄却サンプリング

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

目的とする確率分布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)のサンプルとして採択し、そうでなければ棄却する。

となる。わお!これはあの面積比から円周率を求める計算と同じ!ここから得られるヒストグラムを作ってみると、

mc-rejection

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

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

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

ただし、MCMC ではサンプル間に相関があることには注意が必要。

詳細つり合い条件

遷移確率の決め方として一般的に利用されるのが詳細つり合い条件で、ある確率分布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 国に移住するかを表している。上記の条件が満たされていれば、移住(遷移)の後でも人口が変わらないのがわかると思う。この条件を詳細釣り合い条件と呼ぶ。 https://qiita.com/kaityo256/items/f05f9914eb0ad16afe05#fn2

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

いろいろな MCMC

メトロポリス法

メトロポリス法では、詳細つり合い条件を満たすように以下の遷移確率を用いる。

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\}}

になる。このとき、

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

という手順で状態を次々と更新していく。

標準正規分布のヒストグラムを作ってみると、

mc-metoro

参考


プロフィール画像
tamaosa

釣りと登山が好き。

Privacy PolicyBuilt with Gatsby, © 2019 Tamaosa