やってみた Markov chain Monte Carlo methods, MCMC , gibbs sampling

マルコフ連鎖モンテカルロ法。2変量正規分布からGibbs Samplingする方法を考えてみました。
式を流してもよくわからないので、行間ゼロで理解できるまで細切れにして書いてみます。

Gibbs Sampling

無茶苦茶一般的に書かれている書籍だと以下みたいになっている。
ステップごとに\(\theta=(\theta_1,\theta_2,\cdots,\theta_n)\)の各要素を、
その要素以外の要素の条件付き確率でランダムに発生させる。
\begin{eqnarray}
\theta^0 &=& (\theta_1^0,\theta_2^0,\cdots,\theta_n^0) \\
\theta_i^{t+1} &\sim& p(\theta_i|\theta_1^{t+1},\cdots,\theta_{i-1}^{t+1},\theta_{i+1}^t,\cdots,\theta_n^t)
\end{eqnarray}

2次元空間であれば、ある点を決める際に、
片方の変数\(\theta_1\)を固定して条件付き確率\(P(\theta_2|\theta_1)\)を最大にするパラメタ\(\hat{\theta_2}\)を決める。
その時点で\((\theta_1,\hat{\theta_2})\)が決まる。
次は変数\(\hat{\theta_2}\)を固定して条件付き確率\(P(\theta_1|\hat{\theta_2})\)を最大にするパラメタ\(\hat{\theta_1}\)を決める。

2次元正規分布であれば、片方の確率変数を固定したときの条件付き確率が1変数の正規分布になるから、
これがすごくやりやすい。
Gibbs Samplingは、このように条件付き確率を計算できる必要がある。

2変量正規分布の条件付き確率

2変量正規分布の片方の確率変数を固定すると1変数の正規分布が出てきます。
山の輪切りにして出てきた正規分布の母数を調べたいのですが、導出が無茶苦茶大変そうです。
出てきた正規分布の母数について式変形すると出てくるようです。こちらを参考にさせて頂きました。
この関係式を利用してひたすら式変形していきます。
\begin{eqnarray}
f(x|y) &=& \frac{f(x,y)}{f(x)}
\end{eqnarray}

今、確率変数\(y=y’\)を固定し確率変数\(x\)の正規分布を考えます。
2変量正規分布の分散共分散行列が\(\sum\)であるとします。
\begin{eqnarray}
\sum &=& \begin{pmatrix}
\sigma_{xx} & \sigma_{xy} \\
\sigma_{yx} & \sigma_{yy}
\end{pmatrix}
\end{eqnarray}
出てくる正規分布の平均は以下となります。
確率変数が\(x\)なのに固定した\(y\)が出てくるのがポイント。
\begin{eqnarray}
\mu’ = \bar{x} + \frac{\sigma_{xy}}{\sigma_{yy}} (y’-\bar{y})
\end{eqnarray}
また、出てくる正規分布の分散は以下となります。
\begin{eqnarray}
\sigma’^2 &=& \sigma_{xx} – \frac{\sigma_{xy}}{\sigma_{yy}} \sigma_{yx}
\end{eqnarray}

Python実装例

実際にPythonコードにしてみた図です。2変量正規分布の条件付き確率さえわかってしまえば
あとはコードに落とすだけです。

実行結果