統計学

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

2022年8月15日

ここでは、マルコフ連鎖モンテカルロ法(MCMC法)について説明していきます。

マルコフ連鎖モンテカルロ法(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} ) \)かで決めます。

 
 
 

▼ お問合せはこちら

-統計学
-

© 2024 Yosshi Labo. Powered by AFFINGER5