前回 : Pythonで学ぶ！ベイズ統計第6回 - ma-hi-ro Hobby

Pythonで学ぶ！ベイズ統計第7回

前回は、PyMC3を使って、MCMCを実行し、パラメータ推定を行った
今回は、MCMCの手順を説明して、実際にPythonで実装・実行してみる
ただ、私個人の認識が多大に含まれています

Pythonで学ぶ！ベイズ統計第7回

MCMCの手順

名前の通り、モンテカルロ法（乱数法）とマルコフ連鎖と用いて、定常分布への収束を行っている
- モンテカルロ法については、円周率を求める方法が分かりやすい
- マルコフ連鎖については、「前日の天気から翌日の天気を予測する例」が分かりやすい
それぞれ興味があれば検索してみてください

実際の手順は下記のようになっている

前提
- 推定対象のパラメータが $\theta$
- $\theta$ についての事前分布が $P(\theta)$
- $data$ が得られた場合の、尤度分布が $P(data | \theta)$
- 事後分布（ $data$ が得られた場合の、 $\theta$ の更新分布） $P(\theta | data)$ を求める

実際の手順は下記
1. まずは、パラメータの適当な初期値 $\theta_0$ を定める
2. パラメータ $\theta_0$ を用いて、カーネルを求める
- $prior = P(\theta_0)$
- $likelihood = P(data | \theta_0)$
- カーネル $posterior \propto likelihood \times prior$
3. 次に、提案分布から、乱数を使って、提案値 $\theta_1$ を求める
4. 提案値 $\theta_1$ での、カーネルも求める
5. ２つのカーネルの比と、「1」から提案値の採択率を求める
6. 採択率から、乱数を使って、採用するかどうかを決める
7. 採用する場合は、 $\theta_0 = \theta_1$ として、3. ~ 6を繰り返す
8. 採用しない場合は、 $\theta_0$ のままで、再度3. ~ 6を繰り返す

実際は、提案分布や、採択率の計算方法などで、効率や精度を高められるらしい

交通事故の例1

実際に前回の交差点での交通事故の例を実装版MCMCで求めてみる

まずは、PyMC3モジュールのインポート

import numpy as np
from scipy import stats

data_s = [4]

alpha = 1.6
beta = 1.0

np.random.seed(1)

lambda0  = np.random.gamma(shape=alpha, scale=1.0/beta)
prior0   = stats.gamma.pdf(x=lambda0, a=alpha, scale=1.0/beta)
likelihood0 = np.prod(stats.poisson.pmf(data_s, lambda0))
posterior0 = likelihood0 * prior0

lambda_s = [lambda0]
for i in range(10000):
    while True:
        lambda1 = np.random.normal(lambda0, 1.0)
        if lambda1 >= 0.0:
            break

    prior1   = stats.gamma.pdf(x=lambda1, a=alpha, scale=1.0/beta)
    likelihood1 = np.prod(stats.poisson.pmf(data_s, lambda1))
    posterior1 = likelihood1 * prior1

    p_move = min(1, posterior1 / posterior0)
    uni_random = np.random.uniform()
    if p_move >= uni_random:
        lambda0    = lambda1
        posterior0 = posterior1
    else:
        lambda0    = lambda0
        posterior0 = posterior0
    lambda_s.append(lambda0)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.hist(lambda_s[2000:], density=True, bins=20, range=[-0.25, 9.75], rwidth=0.6)
ax2.plot(lambda_s, '-')

ax1.grid(); ax1.set_ylabel(r'$P$'); ax2.set_xlabel(r'$\lambda$')
ax2.grid(); ax2.set_ylabel(r'$\lambda$'); ax1.set_xlabel(r'$n$')

lambda_s = np.array(lambda_s)
print('%7.3f %7.3f' %(lambda_s[2000:].mean(), lambda_s[2000:].std()))

  2.847   1.163

ほぼほぼ前回と同じ結果

交通事故の例2

データ数を増やした場合

import numpy as np
from scipy import stats

data_s = [4, 3, 1, 2, 1, 5]

alpha = 1.6
beta = 1.0

np.random.seed(1)

lambda0  = np.random.gamma(shape=alpha, scale=1.0/beta)
prior0   = stats.gamma.pdf(x=lambda0, a=alpha, scale=1.0/beta)
likelihood0 = np.prod(stats.poisson.pmf(data_s, lambda0))
posterior0 = likelihood0 * prior0

lambda_s = [lambda0]
for i in range(10000):
    while True:
        lambda1 = np.random.normal(lambda0, 1.0)
        if lambda1 >= 0.0:
            break

    prior1   = stats.gamma.pdf(x=lambda1, a=alpha, scale=1.0/beta)
    likelihood1 = np.prod(stats.poisson.pmf(data_s, lambda1))
    posterior1 = likelihood1 * prior1

    p_move = min(1, posterior1 / posterior0)
    uni_random = np.random.uniform()
    if p_move >= uni_random:
        lambda0    = lambda1
        posterior0 = posterior1
    else:
        lambda0    = lambda0
        posterior0 = posterior0
    lambda_s.append(lambda0)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.hist(lambda_s[2000:], density=True, bins=20, range=[-0.25, 9.75], rwidth=0.6)
ax2.plot(lambda_s, '-')

ax1.grid(); ax1.set_ylabel(r'$P$'); ax2.set_xlabel(r'$\lambda$')
ax2.grid(); ax2.set_ylabel(r'$\lambda$'); ax1.set_xlabel(r'$n$')

lambda_s = np.array(lambda_s)
print('%7.3f %7.3f' %(lambda_s[2000:].mean(), lambda_s[2000:].std()))

  2.507   0.591

こちらもほぼほぼ前回と同じ

初期値をずらした場合

あえて初期値を20から始めた場合

import numpy as np
from scipy import stats

data_s = [4, 3, 1, 2, 1, 5]

alpha = 1.6
beta = 1.0

np.random.seed(1)

lambda0  = 20.0
prior0   = stats.gamma.pdf(x=lambda0, a=alpha, scale=1.0/beta)
likelihood0 = np.prod(stats.poisson.pmf(data_s, lambda0))
posterior0 = likelihood0 * prior0

lambda_s = [lambda0]
for i in range(10000):
    while True:
        lambda1 = np.random.normal(lambda0, 1.0)
        if lambda1 >= 0.0:
            break

    prior1   = stats.gamma.pdf(x=lambda1, a=alpha, scale=1.0/beta)
    likelihood1 = np.prod(stats.poisson.pmf(data_s, lambda1))
    posterior1 = likelihood1 * prior1

    p_move = min(1, posterior1 / posterior0)
    uni_random = np.random.uniform()
    if p_move >= uni_random:
        lambda0    = lambda1
        posterior0 = posterior1
    else:
        lambda0    = lambda0
        posterior0 = posterior0
    lambda_s.append(lambda0)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(9, 4))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)

ax1.hist(lambda_s[2000:], density=True, bins=20, range=[-0.25, 9.75], rwidth=0.6)
ax2.plot(lambda_s, '-')

ax1.grid(); ax1.set_ylabel(r'$P$'); ax1.set_xlabel(r'$\lambda$')
ax2.grid(); ax2.set_ylabel(r'$\lambda$'); ax2.set_xlabel(r'$n$')
ax2.set_xlim(-10, 200)

lambda_s = np.array(lambda_s)
print('%7.3f %7.3f' %(lambda_s[2000:].mean(), lambda_s[2000:].std()))

  2.518   0.591

100回ほどで収束している

おまけ計算で求める

結局、事後分布の確率（カーネル）を比較して、パラメータの事後分布を得ているのなら、
直接事後分布を計算して、規格化すれば良いのでは？と思い、実際に計算してみた

$\lambda$ を 0.1~10まで振って、そのときの事前確率と尤度から事後分布（相対値）を求める

import numpy as np
from scipy import stats

data_s = [4, 3, 1, 2, 1, 5]

alpha = 1.6
beta = 1.0

lambda_calc_s = []
P_s = []
for lambda_ in np.arange(0.1, 10.0, 0.01):
    prior1   = stats.gamma.pdf(x=lambda_, a=alpha, scale=1.0/beta)
    likelihood1 = np.prod(stats.poisson.pmf(data_s, lambda_))
    posterior1 = likelihood1 * prior1

    lambda_calc_s.append(lambda_)
    P_s.append(posterior1)

P_s = np.array(P_s)
P_s = P_s / np.sum(P_s)

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(5, 4))
ax1 = fig.add_subplot(1, 1, 1)

ax1.plot(lambda_calc_s, P_s, '-')

ax1.grid(); ax1.set_ylabel(r'$P$'); ax1.set_xlabel(r'$\lambda$')