Gibbs samplingのコードを書いてみた #Python

ほぼ同等の内容（そして説明も簡潔！）を記載している【M-H法】Python で MCMC を書いてみた【Gibbs-Sampler】があるにも関わらず、Gibbs samplingのサンプルコードを書いてみました。
（限りなく、車輪の再発明感満載ですけども（涙））

何故公開？

冒頭に公開されているpythonコードがありますが、何故公開？
兎にも角にも話がややこしいMCMCなので、初学者向けに「取り敢えず動くコード」があると考えの足しになるかなー、と恩着せがましい思ったことがキッカケです。

ですので、「書いてみました」というよりかは、自分で書いたコードの量が意外に短く済んだので、「何だかよく分からんが、結局何なんだ！」というフラストレーションをお持ちの方の少しでも足しになれば、、、

MCMCとは？

こちらも先人のMCMCについて整理してみた。からどうぞ。

今回生成する分布

2次元正規分布ですねー：
$$
p(x,y)=\frac{1}{Z}\exp\big(-\frac{x^2-2bxy+y^2}{2}\big),
Z=\frac{2\pi}{\sqrt{1-b^2}}
$$
下記のpythonコード内で使用する変数に合わせました。$Z$は正規化定数。

pythonコード

## 準備
import seaborn as sns
import numpy as np
from scipy.stats import norm

## 2次元正規分布の相関を決める定数
b = 0.8
## サンプルサイズの決定
N = 10000

## sampling出来たものを入れとくリスト
x = []
y = []

## 初期点はこんな感じ
x_init = 3.0
y_init = 9.0

## 2次元正規分布に従う乱数発生。
## for文のelse以降が、俗に言う「カクカク」動くところ
for i in range(N):
    if i==0:
        x.append(x_init) ##初期点(x座標)
        y.append(y_init) ##初期点(y座標)
    else:
        x.append(norm(loc=y[i-1]*b, scale=1.0).rvs(size=1.0)[0]) #y座標を固定し、x座標をN(中心点=固定したy座標,1)から選択
        y.append(norm(loc=x[i]*b, scale=1.0).rvs(size=1.0)[0]) #x座標を固定し、x座標をN(中心点=固定したx座標,1)から選択

%matplotlib inline ##Jupyter Notebookで描画したい場合はこのオマジナイ！
sns.jointplot(np.array(x),np.array(y))

初期値関連でひとこと

初期値を変えても、そのうち定常分布になるよー、という文脈なので、「定常分布」だけ欲しい場合は初期値周辺の少し経った辺りのサンプリングを切ってしまうことも考慮ポイントですが、
今回は「兎も角、シンプルなコードを！」をモットーにしましたw

Reference

MCMC講義（伊庭幸人）難易度★★：MCMCを本当に平易なところから解説されている動画です。
【統計学】マルコフ連鎖モンテカルロ法(MCMC)によるサンプリングをアニメーションで解説してみる。:同じ分布を丁寧にアニメーションまで付加してらっしゃいます。