Gibbs sampling 코드를 써봤어요.

기본적으로 같은 내용(그리고 설명도 간단합니다!)[M-H법] 파이톤으로 MCC를 써봤어요. [Gibs-SApler]의 기록에도 불구하고Gibssampling의 샘플 코드를 써 보았습니다.
(무한, 바퀴의 재발명감 가득(눈물))

왜 공개해요?

처음에 공개된python 코드가 있는데 왜 공개합니까?
아무래도 MCMC가 복잡하기 때문에 초보자들에게'잠시 이동하는 코드'가 있다는 생각은 고민을 돕는 것이 아닐까.
그래서'써봤어'라기보다는 본인이 쓴 코드의 양이 의외로 짧아져서 "어떻게 된 일인지 결과가 뭐였는지!"이런 좌절감을 느끼는 사람이 조금이라도 있다면

MC란 무엇인가?

이것도 선인의 MC 정리 좀 했어요.부터다.

이번에 생성된 분포

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$는 귀일화 상수입니다.

개인 코드

## 準備
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))

초기값과 관련된 말

초기 값이 바뀌어도 머지않아 안정적 분포로 바뀔 수 있다. 따라서'안정적 분포'만 원한다면 초기 값 근처에서 시간이 조금 지나면 샘플링을 차단해야 한다는 것을 감안하면
이번에는 "토끼도 뿔, 간단한 코드!"구호

Reference

MCMC과정(이정행인)난이도★: MC를 진정으로 간단명료하게 설명하는 영상입니다.

[통계학] 애니메이션으로 마코프 체인 몬테카로법(MCMC)의 샘플링을 해석한다.: 동일한 분포를 애니메이션에 세밀하게 추가합니다.