機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングス法(2)

今回は多変量正規分布をMCMC法でサンプリングしてみました。先回書いたコードを少し修正するだけなので実装は簡単です。多変量正規分布で書いた分散共分散行列が


\Sigma 
=
 \left(
    \begin{array}{cc}
      1 & 0.7 \\
      0.7 & 1
    \end{array}
  \right) \tag{1}

のものをサンプリングしてみます。平均は0としました。また、提案分布も先回同様に正規分布を用いたものにしました。

f:id:opabinia2:20180410211230p:plain うまくサンプリングできていますね。先回同様、提案分布の分散が適切でないと、以下のようにサンプリングがうまくいかない場合があります。分散が小さすぎて分布全体からサンプリングできていません。

f:id:opabinia2:20180410211237p:plain

今回のコードです。多変量正規分布だと面白くないのでもうちょっと違う分布も試してみたいですね。面白い形の関数ないかな。

# マルコフ連鎖モンテカルロ法(メトロポリスヘイスティングス法)のテスト
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm


# 目標分布(多変量正規分布)の分散共分散行列
s = np.array([[1, 0.7], [0.7, 1]])


# 提案分布
def q(x):
    mu= [0, 0]
    cov = [[0.5, 0], [0, 0.5]] 
    z = np.random.multivariate_normal(mu, cov).T
    return x + z


# サンプルしたい関数
def p(x1, x2):
    x = np.array([x1, x2])
    return (1/np.sqrt((2*np.pi)**2*np.linalg.det(s))) * np.exp(-0.5*np.dot(np.dot(x.T,np.linalg.inv(s)),x))


# MCMC
def MCMC(x):
    while True:
        # サンプル候補
        z = q(x)
        # 受理するまでサンプル候補をとる
        if p(z[0], z[1])/p(x[0], x[1]) > np.random.uniform(0, 1):
            return z


vec_p = np.vectorize(p)

# ランダムシードを固定
np.random.seed(0)

# サンプル数
N = 3000

# バーンイン
B = 1000

# 初期値
x = np.array([2, 2])

samples = np.empty([N, 2])

# MCMC
for i in range(N):
    x = MCMC(x)
    samples[i, :] = x


x, y = np.meshgrid(np.linspace(-3, 3, 100), np.linspace(-3, 3, 100))
z = vec_p(x, y)
plt.contourf(x, y, z, 100, cmap=cm.coolwarm)
plt.scatter(samples[B:, 0], samples[B:, 1],c="green", marker="*",alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.show()