統計学

わかりやすいモンテカルロ法とマルコフ連鎖モンテカルロ法(MCMC法)

2022年8月15日


モンテカルロ法とは、
シミュレーションや数値計算を乱数を用いて行う方法の総称で、総称なのでさまざまなモンテカルロ法が世の中には存在しています。
さまざまな計算を乱数を用いて処理するような仕組みを提供するものになります。
 

よくある方法としては、円周率を計算する方法。
さらには棄却・受容サンプリングを用いると、サンプリングを用いて不明な分布のサンプリングを行います。
さらには重点サンプリング法を使うと、ある解析的に計算ができない積分計算を乱数を用いて近似計算を行うことができます。(これをモンテカルロ積分という)

さまざまな計算や課題などをサンプリングや乱数を用いて解決する手法全般のことですね!
 

モンテカルロ法

モンテカルロ法は、
円を囲む正方形を用意する。
この正方形にランダムに乱数を出していく。

そして、とにかく乱数を乱れ打ちしていく。隙間無くなるくらいに
もちろん重なりもなってよし!
で正方形の中に円があって、その円の面積を求めたいとする。
そこで乱れ打ちの結果できた点を数えます。
 

正方形全体は乱数発生させた回数であり、\(N\)回
そしてその円の中に撃たれた回数が\(X\)回とすると、
正方形と円の撃たれた回数の比は\(X:N\)
つまり、これは正方形と円の面積比でもある。\(X:N\)
つまり、面積と当たった球の数は比例するって話ですね。
面積が大きければその分当たりやすくなるし、小さければその分当たりづらくなるしということです。
 

よって正方形の面積を\(S\)とすれば、円の面積は\(S \cdot \displaystyle \frac{X}{N}\)になる。
よって正方形の半径が\(R\)とすれば、逆算して\(Rπ = S \cdot \displaystyle \frac{X}{N}\) により、\(\pi\)がもとまる。
 

大数の法則によって、収束していくので、
数が多ければ精緻になっていく。
円の中隙間なく全体に撃たれるまではやる必要はなく、
正直少なくても問題はないが、その場合運良くそこに当たって、
円の中に100回、正方形と円の間の部分に3回とかだと面積比が100:103になるので明らかおかしい。
なので、大数の法則から数を増やせば増やすほどその運よくが調整されて、精緻になっていきます。
 

実際にプログラムで発生させてみる

こうなると難しい面積なども容易に求めることができそうである。

以下の方法はまさにこの乱数発生して、面積でサンプリングをしていく。
これが基本となる。
 

棄却サンプリング

ここでは重点サンプリングについて説明しますが、
マルコフ積分と呼ばれる解析的に計算をすることのできない積分をサンプリングを用いて推定するという仕組みを提供します。
 

pythonや統計解析ソフトRなどでは正規分布などよくある基本的な関数からのサンプリングはできますが、
不明な関数からのサンプリングはできない場合に、どのようにしてサンプリングを行うかが争点となります。
 

棄却サンプリングとは,サンプルを生成したい分布を覆うようなよく知られた分布から生成したサンプルを用いて,サンプリングを行う方法です.
予測分布\(p(z)\)のサンプルを提案分布\(q(z)\)から求めることを考えます.
このとき\(p(z) < kq(z)\)が満たされているものとします
→ つまり、\(p(z) < kq(z)\)なので、\(kq(z)\)でサンプリングしたい\(p(z)\)を覆い尽くすようにする。
まさにこれって、モンテカルロ法の正方形で囲むものと同じ。
 

まず知っている\(q(z)\)からサンプリングをします。
そして\(p(z)\)と\(q(z)\)の比率は0から1なので、その間に対応して、0から1でランダムな数値を出します。
そして\(q(z)\)からもデータのサンプリングをします。
この2つの結果を用いて以下の判定をします。
\(u < \displaystyle \frac{p(z)}{cq(z)}\)であれば、
 

ただ注意としては、なかなかこの\(p(z) < kq(z)\)を満たす、全ての\(z\)に対してこの式が成り立つような\(q(z)\)を満たす関数を出すのが難しい。 なので、この棄却サンプリングに加えてさらに条件をつけて、サンプリングをするような手法が次のやり方になります。  
 

でも今回結果的には\(\displaystyle f \left(x | u < \frac{p(z)}{cq(z)} \right)\)からサンプリングして上記を満たすサンプリングをしたわけだが、本当にこれだけで\(p(z)\)からのサンプリングに帰着できるのでしょうか?

まずベイズの定理から以下のように変換することができます。\(\theta\)はパラメータ群を表しているとして、
\begin{eqnarray}
f(\theta|x)
&=& \displaystyle \frac{f(x,\theta)}{f(x)} \\
&=& \displaystyle \frac{f(x|\theta)f(\theta)}{f(x)} ・・・① \\
\end{eqnarray}
ここで\(f(x)=\displaystyle  \int f(x|\theta)f(\theta)d\theta \)であることを利用して、
\begin{eqnarray}
f(\theta|x)
&=& ①\\
&=& \displaystyle \frac{f(x|\theta)f(\theta)}{\displaystyle \int f(x|\theta)f(\theta)d\theta }
\end{eqnarray}
となります。

これを利用して、以下のように変換できます。

実際には以下のように帰着でき、
\begin{eqnarray}
f\left(z|u < \displaystyle \frac{p(z)}{cq(z)} \right)
&=& \frac{P\left(u < \displaystyle \frac{p(z)}{cq(z)}|z\right)q(z)}{\displaystyle \int P\left(u < \displaystyle \frac{p(z)}{cq(z)}|z\right)q(z)dz} \\
&=& \frac{\displaystyle \{p(z)/cq(z)\}q(z) }{\displaystyle \int \{p(z)/cq(z)\}q(z)dz} \\
&=& \displaystyle \frac{ p(z) }{ \displaystyle \int p(z)dz} \\
&=& p(z)
\end{eqnarray}
上の式の意味としては、条件付き確率なので、
まずサンプリングをした上で、\(u\)より大きいかの確率を出してます。
その中で計算していくと、そのサンプリングは結果的には\(p(z)\)から出たものと同値なので、このサンプリング方法でうまくいくということになります。

 

ただ欠点としては、そういった\(q(z)\)を求めることが難しいです。
もし\(q(z)\)を大きく作るのでも良さげ!って思いますが、その場合、各\(z\)ごとの\(p(z)\)と\(q(z)\)の比率がめちゃくちゃ小さくなるので、
0から1のサンプリングをするとき、受容するデータが極端に少なくなり、そうなると処理が長時間に及ぶ可能性があるためです。
 

重点サンプリング

確率分布を求めることは難しいが、
確率変数に対して確率を求めることはできる場合、
あえて、期待値で考えることで求めれる。
 

\begin{eqnarray}
H = \int h(x)dx
\end{eqnarray}
が解析的に求めれない場合があるとする。この計算を求めるためのサンプリングになります。
ここでh(x)を覆うような新たな関数\(\pi(x)\)とします。\(\pi(x) は0を取らない\)とした時、
( \(\pi(x) > h(x)\) )
\begin{eqnarray}
H
&=& \displaystyle \int \frac{h(x)}{\pi(x)}\pi(x)dx \\
&=& \displaystyle E_{\pi} \left[\frac{h(x)}{\pi(x)}\right]・・・①
\end{eqnarray}
上の2行目を見ると、期待値の定義の形式に似ています。
通常\(E[X] = \int xf(x)dx\)です。
なので、期待値の定義の意味は、f(x)から生成された確率変数\(x\)の期待値という意味になるので、
上の\(①\)は、\(\pi(x)\)から生成された確率変数\(\displaystyle \frac{h(x)}{\pi(x)}\)の期待値ということになります。
こうすることで、解析的に求めれなかったh(x)は確率変数の方に持っていくことができたのが分かります。
 

なぜこうする必要があり、こうすることによるメリットは何か?
それは平均と期待値の意味するところの違いになります。
期待値とは「統計学で、確率の見地から算定した平均値」で、
平均は「各値の総和を個数で割った値」のことです。
つまり、期待値は確率を用いて計算をした平均のことなので、
今上記積分の期待値を計算できない今、あと残りの平均を求める方法としては実際に値の総和を出して、それを個数で割って計算をするということに帰着できそうです。

 

そして、確率変数の和は中心極限定理と大数の法則によって収束することが知られています。
その原理を用いることで、
\begin{eqnarray}
H &=& \displaystyle E_{\pi} \left[\frac{h(x)}{\pi(x)}\right]・・・① \\
&\rightarrow& \frac{1}{n} \sum_{i=1}^{n} \frac{h(x_{i})}{\pi(x_{i})}
\end{eqnarray}
元々の仮定としては、\(H\)を解析的に計算はできない状態で、\(h(x)\)はわかっている状態です。
さらに\(\pi(x)\)も自分で選択した関数なのでわかっているはずなので、
\(h(x)\)と\(\pi(x)\)のxにランダムで値を入れていき、上の平均の計算をしていくことで、最終的なHを求めることができるということになります。
上記の平均式を見ると、分母は\(\pi(x)\)なので、この値が極端に大きすぎる、つまり覆いかぶしても極端に\(h(x)\)から離れていては計算が計算が小さくなってしまいます。

 

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

ここまでで解説したマルコフ連鎖とモンテカルロ法を組み合わせて、事前分布と尤度分布から、事後分布をシミュレートするのがマルコフ連鎖モンテカルロ法です。

マルコフ連鎖モンテカルロ法は、ベイズの定理の事後分布を導出するためのもの。
マルコフ連鎖モンテカルロ法はそれ自体サンプリングではなく、事後分布を求めるサンプリングの総称。
このMCMC法に分類されるものとしては、

MCMCを使う目的は、統計モデルのパラメタを推定することです。

まずモンテカルロ法は上で説明済み。
マルコフ連鎖は、1つ前の結果を用いて、次に出る値の確率を予測するものになります。となると、上記のモンテカルロ法では1つ前の結果などは関係せずあくまで毎回の試行結果のみを判定するが、マルコフ連鎖は1つ前の結果を揉みつつサンプリングをしていくところに違いがあります。
マルコフ連鎖は以下で扱ってますので、ぜひ。
マルコフ連鎖

モンテカルロ法でデータのサンプリングをして、マルコフ連鎖で1つ前のサンプリング結果から今回の結果の確率を用いたハイブリット方式ということになります。

マルコフ連鎖は、1つ前の状態から次の状態を予測するためのもの。
モンテカルロ法は、上で説明した通り、サンプリングや乱数で課題を解決する。

マルコフ連鎖モンテカルロ法はなのでこの2つを混ぜた手法になります。
1つ前の状態から次の状態を予測する際に、サンプリングや乱数などを用いて計算するというイメージを持ってもらえればと思います。
 

主に、有名な手法として、

  • ギブスサンプリング
  • メトロポリス・ヘイスティング法

の2つがあります。

となるとこの2つを合体させたものがいわゆるMCMC法になります。

通常ベイズの定理で、
事後分布 = 尤度関数×事前分布
から解析的に事後分布を計算して求めることができない場合があります。
共役事前分布や共役事後分布が同じ分布であったりすると、簡単に求めることができたりはしますが。
そしてパラメータを求めるmode(最頻値)を求めるが、その場合事後分布であればそのまま盛り上がってるところで最頻値パラメータになります。
しかし解析的に事後分布がわからないのであれば、サンプリングをして最も頻度が多く出たパラメータを採用するということになります。
そもそも事後分布わからなくても次の事前分布で使用するパラメータはわかるから!
 

マルコフ連鎖
\begin{eqnarray}
K(x_{n+1}|x_{n},x_{n-1},...,x_{0}) = P(x_{n+1}|x_{n})
\end{eqnarray}

ギブスサンプリング

モンテカルロ法はもう何でも受け入れるという形だったが、
ギブスサンプリングはある一定の基準を設けて、その円の中に点を打つかどうかを決める。
これは、事後分布が決まっている状態で、その事後分布からパラメータを解析的に求めれない場合に使用するサンプリング法。

でも事後分布も解析的に求めれない場合がある。そうなったらサンプリングできないよね。
その対策としてMH法がある

 

マルコフ連鎖では以下のような性質がある。

マルコフ連鎖は正則条件の下では、反復することによって確率標本の分布が不変分布に収束するという性質があるので、
ベイズ更新の事後分布に対してマルコフ連鎖を都度行っていくことによって、そこで得られた標本が不変分布に収束する、つまりこれは求めているものになります。

 
パラメータを求める必要があるわけで、
例えば、パラメータが\(\mu, \sigma^{2}\)だとして、確率変数は\(x\)とした時、
それぞれ交互にサンプリングしていく方法がギブスサンプリングになります。
当然この2つのパラメータが同時分布\(P(\mu, \sigma^{2}|x)\)からはサンプリングはできません。
なので、それぞれ\(\mu\)のみの関数、\(\sigma^{2}\)のみの関数にしていきます。
要は周辺分布ということですね!
なので、
\(P(\mu|x, \sigma^{2})\)
\(P(\sigma^{2}|x, \mu)\)
を求めます。
\begin{eqnarray}
\mu^{1} → (\sigma^{2})^{1} → \mu^{2} → (\sigma^{2})^{2} → ...
\end{eqnarray}
 

予測分布

もちろん、MCMC法の1つであるこのギブスサンプリングは事後分布はパラメータによる関数なので、このサンプリングを用いることでパラメータを得ることができます。
そして事前分布にて、その得たパラメータの出る確率を求め、そして尤度関数に適用します。
 

予測分布を定義します。
今得られたデータ\(x_{n}\)から、次のデータ\(x_{n+1}\)を出したいとします。
この時、パラメータは\(\mu, \sigma^{2}\)とすると、
予測分布は以下のように設定できます。
\begin{eqnarray}
f(x_{n+1}|x_{n}) &=& \int\int f(x_{n+1}|\mu, \sigma^{2})f(\mu, \sigma^{2}|x_{n})d\mu d\sigma^{2}
\end{eqnarray}
一番左はそれこそ次の\(x\)を求める予測分布。一番右は今得られているxから求めた\(\mu\)と\(\sigma^{2}\)なので事後分布。
\(x\)→\(\mu\)と\(\sigma^{2}\)→次の\(x\)を求めるという形です。
そして2つのパラメータを積分させることで、この2つのパラメータがとりうる全ての値で計算をして合算した値がこの予測分布になります。
実際この2つのパラメータで積分すれば右の項の真ん中は消えるので、
\(f(x_{n+1}|\mu, \sigma^{2})f(\mu, \sigma^{2}|x_{n})\)は、間の\(\mu, \sigma^{2}\)は消えて、\(f(x_{n+1}|x_{n})\)になります。
 

予測分布の求め方\(①\)

実際これを計算して求めようとするとまずほぼほぼ無理です。
なのでどうするかというと、上記の積分の意味を考えると全ての\(\mu\)と\(\sigma^{2}\)の組み合わせでの和なので、
上のギブスサンプリングで得られた\(\mu\)と\(\sigma^{2}\)の組み合わせを以下のように定義します。
\((\mu^{(i)}, {\sigma^{2}}^{(i)})\)
iはサンプリングで得られた個数です。5個得たのであれば、\(i=1,2,3,4,5\)で、5組得られたことになります。
 
よって予測分布は以下のように表せます。
\begin{eqnarray}
f(x_{n+1}|x_{n}) &=& \sum_{i=1}^{5} f(x_{n+1}|\mu^{(i)}, {\sigma^{2}}^{(i)})f(\mu^{(i)}, {\sigma^{2}}^{(i)} |x_{n})
\end{eqnarray}
 

予測分布の求め方\(②\)

もっと簡単に考えると、
そもそも\(x_{n}\)から\(\mu\)と\(\sigma^{2}\)は求まっています。
なので、であれば、この\(\mu\)と\(\sigma^{2}\)の組み合わせから、(上の組み合わせと同じ意味)
直接、\(x_{n+1}^{(i)} \sim N(\mu^{(i)}, {\sigma^{2}}^{(i)})\)なので、
直接それぞれの\(i\)の組み合わせで、サンプリングをして、
ヒストグラムを書いて、一番出る最頻値(mode)を採用すればいい。

なので、予測分布の\(x_{n}\)を条件としたことは、間接的に\(x_{n}\)を用いて\(\mu\)と\(\sigma^{2}\)は求まっているので、
サンプリングでヒストグラム書いて求めればそれでいい!
 

またこの予測分布の関数が正規分布などよく知られているような確率分布でないのであれば、サンプリングは少々難しいです。
有名な分布であれば、pythonやRで簡単にサンプリングができます。
知られていない場合は、それこそ上で説明した、モンテカルロ法を用いてサンプリングを行います。
「モンテカルロ法はさまざまな確率分布からの乱数を発生させる手法」なので。
 

メトロポリス・ヘイスティングス法

ある分布の面積を求めたいとする。ただこの分布がわからず、ある分布はわかってる時、

 

あれか、右に行けば行くほど確率は低くなるからか、、
だから右側から真ん中への遷移は確率高いけど、
真ん中から右への遷移は確率が低いので、自分のところにとどまるみたいな話か。。。
でそれを採択するかどうかは、右にいけばいくほど採択はなるべくしたくないので、
高確率で真ん中にとどまろうとする。でも少しの確率で右に移動するかもしれないってことか。
 

つまり一様分布であればどこも一緒だけど、
正規分布とかだとその枠周り、ギリギリ周りだとほんと確率が小さいのでそこでは採択はできる限りしたくないみたいな話かー?

ここでマルコフ連鎖を使うという考えのもと、今状態\(x\)から\(x'\)へ遷移する際の遷移核を\(K(x,x')\)とします。
そうすると、マルコフ連鎖と同じようなイメージで、\(p(x') = K(x',x)p(x)\)となります。
この遷移核を定義した上で以下の説明をしていきます。
 

マルコフ連鎖の性質「詳細釣り合い条件」

マルコフ連鎖は、そもそも状態が離散の場合の遷移行列になります。なので連続の場合はまた別に定義が必要になります。
ここで連続の場合の遷移行列を\(x\)から\(x'\)に状態が遷移するとして、\(K(x,x')\)とします。
この時、以下の式が成り立つ場合、マルコフ連鎖が不変分布をもち、収束するということが言えます。

\begin{eqnarray}
K(x,x')p(x') = K(x',x)p(x)
\end{eqnarray}

上記の式の意味はなんとなくわかると思います。
左辺であれば、現在状態\(x'\)にいるとしたとき、別の状態\(x\)に遷移する時の確率を意味していて、
右辺であれば、現在状態\(x\)にいるとしたとき、別の状態\(x'\)に遷移する時の確率を意味しています。
となると、その遷移した結果の確率が一緒であれば、釣り合うのはなんとなく想像できると思います。
 

さらにここでもう少し理解するために、
つまり分布の端にいけば行くほど、そもそもの確率が低くなる分、遷移行列の確率は釣り合うために相対的に高くなるので、押し戻される形になります。(遷移核の確率が高くなるということは押し戻す確率が高くなる)
 

となると、想像できそうなのが、
以下2つです。

  • \(①\) \(K(x,x')p(x') = K(x',x)p(x)\)の式をみると、尤度関数と事前分布の形に似てること
  • \(②\) 分布の端に行くほど、押し戻される可能性が高くなるのであれば、押し戻されるので、これが求めたい事後分布のサンプリングに使えるのでは?さらに確率が高い\(p(x)\)であれば、遷移核の確率は低くなるのでその場を移動しなくなる可能性があるので、そうなるとよく出る最頻値にサンプリングが多く集まりそう)

 

まず右にいる場合、分布の端に行くほど、その状態が出る確率p(x)が小さくなる。そうなると均衡を保つために、自ずと遷移核が大きくなり、戻されやすくなる。
上図の\(x'\)の状態の出る確率p(x)が高くなると、均衡を保とうとするために、自ずと遷移核が小さくなり、その状態から移動しずらくなる。その状態あたりに留まろうとする。
そのため、今いるところから状態の出る確率がより高いところにどんどん遷移していきやすくなります。
そうなると、最終的に一番確率の高いところに収束していくのが想像できそうですが、そうなると満遍なくサンプリングができているとは言えません。
 

そのため、今のところから新しい状態の出る確率がより高いところに遷移していくのは前提の現象として、たまには調整をして、無理やり遷移しづらくして満遍なくいろんな状態(確率変数)に遷移していけるようにする仕組みが必要です。
そのためには、
\(\alpha(x,x') = \min \left{ \displaystyle \frac{K(x,x')p(x')}{K(x',x)p(x)} , 1 \right} \)
とおき、この確率によって、次の遷移先の値を採用するかどうかを判定します。
minなので、2つの値のうち小さい方を採用するわけですが、1は定数、そしてもう1つの値は変わります。
分母 > 分子であれば、1より小さくなるのでその値が採用され、アルファの値は1より小さくなり、その確率によって採用するかどうかを決めます。
逆に分母 < 分子であれば、1より大きくなるので1が採用され、アルファの値は1になるので、必ず採用されます。 分母 > 分子の場合は、遷移先が今よりもよく出る確率変数に遷移する時になります。よく出るからこそ、遷移がそこに留まらないように確率(アルファ)が1より小さくしている。
逆に分母 < 分子の場合は、遷移先が今よりも出にくい確率変数に遷移するときになります。滅多に出ないからこそ、サンプリングする際は確率(アルファ)1で採択しようということです。   離散のマルコフ連鎖
マルコフ連鎖では行列\(P\)によって、別の\(x_{t}\)から\(x_{t+1}\)へ遷移させます。まさに\(P\)は推移核ですね!!!\(x_{t+1} = P x_{t} \)
連続のマルコフ連鎖
連続のマルコフ連鎖では、推移核\(K(x_{t+1}, x_{t})\)によって、\(x_{t}\)から\(x_{t+1}\)へ遷移させます。まさに\(K\)は遷移核ですね!
つまり、基本的にはマルコフ連鎖は離散、連続ともに正則条件によって、\(n \rightarrow \infty \)によって不変分布を構築することができます!
これがいわゆる事後分布になりますね!
ただ、メトロヘイスティング法では、このKではなく、アルファで遷移させる\( (x_{t} \rightarrow x_{t+1}) \)か、その場所に止まる\( (x_{t} \rightarrow x_{t} ) \)かで決めます。

 
 
 

▼ お問合せはこちら

-統計学
-

© 2022 Yosshi Blog Powered by AFFINGER5