「基礎からのベイズ統計学」の例を実装しながら理解する（①MH法まで）

「基礎からのベイズ統計学」で紹介された例をPythonで実装してみる。この記事は基本的には備忘録であり、書籍で詳しく解説してある部分は省略し、ポイントと実装して確認した結果のみ記載する。

作者:豊田秀樹
発売日: 2015/06/25
メディア: 単行本

MCMCを使う理由

何はともあれ、ベイズの定理。

$\displaystyle{f(\theta|x)=\frac {f(x|\theta)f(\theta)} {\int f(x|\theta)f(\theta)d\theta}}$

事後分布の計算には分母の積分が必要になるが、解析的に解けないことがほとんど。
そこで、分子の（事前分布）x （尤度）の部分（本書でいう事後分布のカーネル）を利用して、母数の事後分布に従う乱数を発生させることで、母数に関する推論を行うことを考える。そのための手法がMCMC（マルコフ連鎖モンテカルロ法）である。

詳細釣り合い条件

$\displaystyle{f(\theta|\theta')f(\theta')=f(\theta'|\theta) f(\theta)}$

MCMCが定常分布に収束するための十分条件が詳細釣り合い条件である。これが、あらゆる2点間に成り立っている必要がある。
上式は、 $\theta$ と $\theta'$ の確率の比が1:pならば、 $\theta$ から $\theta'$ への移動は $\theta'$ から $\theta$ への移動のp倍起こりやすい、ということを意味する。

メトロポリス・ヘイスティングス法

MCMCのうち、一番単純な手法。適当な提案分布をマルコフ連鎖の遷移核（ $f(\theta|\theta')$ ）として利用する。詳細釣り合い条件を満たすように遷移核を補正していくことで、事後分布に従うサンプルを得られるようになる。

独立MH法

本書の通り、「波平釣果問題」を題材とする。

波平さんの過去10回の釣果は0匹が6回、1匹が3回、2匹が1回だった。波平さんが次に釣りに行くときに、釣果が0匹となる確率を求めよ。ただし、釣果はポアソン分布に従っているとする。また、事前分布は $\alpha=6, \beta=3$ のガンマ分布とする（平均2、標準偏差0.8程度）。

この場合、ポアソン分布の共役事前分布であるガンマ分布を使っているので、事後分布もガンマ分布となり、解析的に確率を求めることが可能である。ここでは、サンプリング結果の正しさを確認するために、敢えてこの問題を題材としている。

独立MH法では、提案分布を無条件分布とする。つまり、提案される候補値は互いに独立となる。
本書では平均1.0 、分散0.5の正規分布を提案分布としている。この場合、負の値がサンプリングされることがあるが、ポアソン分布の母数は正の値である必要があるため、不適である。今回の実装では提案された値が負の場合は、事後分布のカーネルが0となるようにしている。こうすることで、負の値が提案された場合は必ず受容されないことになる。

独立MH法の欠点として、適切な提案分布を選ばないと収束に時間がかかる、または現実的な時間では収束しないことが挙げられている。
提案分布を正規分布として、平均は1.0に固定、標準偏差を0.5, 1.5, 0.1とした場合でサンプリングを試みる。

f:id:kibya:20190331110425p:plain — 目標分布と提案分布

それぞれについて、6000個のサンプリングを行い、最初の1000個をバーンイン期間として破棄して、残りの5000個でヒストグラムを作成した結果を下に示す。
受容率はそれぞれ、53.7%, 20.1%, 39.7%となった。

提案分布の幅が広すぎる場合（標準偏差：1.5）は、受容率が低くなってしまう。また、幅が狭すぎる場合（標準偏差：0.1）は目標分布からの満遍ないサンプルを得ることが難しく、収束させることは難しい。

f:id:kibya:20190331110504p:plain — 独立MH法によるサンプリング結果

ランダムウォークMH法

独立MH法では適切な提案分布の設定が難しいという問題があった。ランダムウォークMH法では候補の提案にランダムウォークを利用することで、この問題を解決している。提案分布に対象な分布を選ぶと、受容確率（補正係数）の式から提案分布の部分が消え、事後カーネルの比となる。

提案分布を $\mu=0$ の正規分布として、 3通りの $\sigma$ でサンプリングを行なった結果が下記。6000個のサンプリングを行い、最初の1000個はバーンイン期間として破棄している。受容率はそれぞれ、49.3%, 92.2%, 14.6%となっている。提案分布の分散 $\sigma_\epsilon$ が適切な値であれば、独立MH法よりも受容率がかなり高くなることが分かる。

f:id:kibya:20190401092145p:plain — ランダムウォークMH法

次に、トレースプロット（各ステップで採択されたパラメータを順番にプロットしていったもの）を見て、それぞれの $\sigma_\epsilon$ の値の場合で何が起こっているかを確認する。なお、提案された値のうち、破棄されたものを赤点で示している。

f:id:kibya:20190401092219p:plain — ランダムウォークMH法　トレースプロット

$\sigma_\epsilon=0.5$ は適切な値であり、少し離れた所にあった初期値から速やかに0.8付近に近づいていることが分かる。一方、 $\sigma_\epsilon=0.05$ と小さすぎる場合は、受容率こそ高いものの、サンプルされる値が少しずつしか変化しないため、収束するには時間がかかる。また、 $\sigma_\epsilon=2.0$ と大きすぎる場合は、受容率が低くなってしまい、やはりサンプリングの効率が悪くなってしまう。

今回のようにパラメータが1つの場合は $\sigma_\epsilon$ の調整は難しくないが、パラメータ数が多くなると $\sigma_\epsilon$ の調整が難しくなり、受容率が低下するという問題がある。
これを解決するための手法が、ハミルトニアンモンテカルロ法（HMC法）である。次回は、こちらを実装してく。

Pythonコードはこちら。Jupyter notebookを前提に書いている。

%matplotlib inline
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import cm
import seaborn as sns
from tqdm import tqdm_notebook
from scipy.stats import norm, poisson, gamma, uniform
sns.set(style="white")

np.random.seed(1)

x = np.linspace(-1, 5, 100)
plt.plot(x, norm.pdf(x, loc=1.0, scale=0.5), label="proposal (mu=1.0, sigma=0.5)")
plt.plot(x, norm.pdf(x, loc=1.0, scale=1.5), label="proposal(mu=1.0, sigma=1.5)")
plt.plot(x, norm.pdf(x, loc=1.0, scale=0.1), label="proposal(mu=1.0, sigma=0.1)")
plt.plot(x, gamma.pdf(x, a=11, scale=1./13), linestyle="--", label="posterior")
plt.xlabel("theta")
plt.ylim(-0.1, 2)
plt.legend()
plt.tight_layout()
plt.savefig("proposal_dists.png")


lam = 2
data = [0, 0, 0, 0, 0, 0, 1, 1, 1, 2]


def posterior_kernel(lam):
    if lam < 0:
        return 0
    likelihood = poisson.pmf(data, mu=lam).prod()
    prior = gamma.pdf(lam, a=6, scale=1./3)
    return prior * likelihood


# 独立MH法
def independent_MH(mu_proposal, scale_proposal, ax, theta_init=1.0, n_sample=6000, burn=1000):
    theta_current = theta_init
    posterior = [theta_current]
    count = 0 # count accepted proposals
    
    for i in tqdm_notebook(range(n_sample)):

        proposal_dist = norm(loc=mu_proposal, scale=scale_proposal)
        theta_proposal = proposal_dist.rvs()

        # Lieklihood x prior
        p_current = posterior_kernel(theta_current)
        p_proposal = posterior_kernel(theta_proposal)

        # probability to accept proposal
        p_accept = proposal_dist.pdf(x=theta_current) * p_proposal / \
        (proposal_dist.pdf(x=theta_proposal) * p_current)

        if np.random.rand() < p_accept:
            count += 1
            theta_current = theta_proposal

        posterior.append(theta_current)
    
    g = sns.distplot(posterior[burn:], ax=ax, label="posterior samples")
    g.plot(x, gamma.pdf(x, a=11, scale=1./13), label="theoretical posterior")
    g.legend()
    g.set_title(f"proposal mu: {mu_proposal}, sigma: {scale_proposal}")
    g.set_xlim(0, 2.5)
    print("accept ratio: {:.2%}".format(count / n_sample))


fig, axes = plt.subplots(1, 3, figsize=(14, 4))
independent_MH(1.0, 0.5, axes[0])
independent_MH(1.0, 1.5, axes[1])
independent_MH(1.0, 0.1, axes[2])
plt.savefig("independent_MH.png")


# ランダムウォークMH法
def randomwalk_MH(scale_proposal=0.5, theta_init=4.0, n_sample=6000, burn=1000):
    theta_current = theta_init
    posterior = [theta_current]
    rejected = []
    count = 0 # count accepted proposals
    
    for i in tqdm_notebook(range(n_sample)):

        # 前のサンプルに基づき、次の候補が提案される
        theta_proposal = norm.rvs(loc=theta_current, scale=scale_proposal)

        # Lieklihood x prior
        p_current = posterior_kernel(theta_current)
        p_proposal = posterior_kernel(theta_proposal)

        # probability to accept proposal
        p_accept = p_proposal / p_current

        if np.random.rand() < p_accept:
            count += 1
            theta_current = theta_proposal
        else:
            rejected.append((i, theta_proposal))

        posterior.append(theta_current)
        
    print("accept ratio: {:.2%}".format(count / n_sample))
    return np.array(posterior), np.array(rejected)


def plot_posterior(posterior, burn=1000):
    g = sns.distplot(posterior[burn:], label="posterior samples")
    g.plot(x, gamma.pdf(x, a=11, scale=1./13), label="theoretical posterior")
    g.legend()
    g.set_xlim(0, 2.5)
    

def plot_trace(posterior, rejected):
    plt.plot(posterior, label="trace")
    plt.scatter(rejected[:, 0], rejected[:, 1], marker="o", color="r", s=4, label="rejected")
    plt.legend()


sigma1, sigma2, sigma3 = 0.5, 0.05, 2.0
posterior1, rejected1 = randomwalk_MH(sigma1)
posterior2, rejected2 = randomwalk_MH(sigma2)
posterior3, rejected3 = randomwalk_MH(sigma3)

fig, axes = plt.subplots(figsize=(14, 4))
plt.subplot(131)
plot_posterior(posterior1)
plt.title(f"proposal sigma: {sigma1}")

plt.subplot(132)
plot_posterior(posterior2)
plt.title(f"proposal sigma: {sigma2}")

plt.subplot(133)
plot_posterior(posterior3)
plt.title(f"proposal sigma: {sigma3}")

plt.savefig("randomwalk_MH.png")


fig, axes = plt.subplots(figsize=(10, 8))
plt.subplot(311)
plot_trace(posterior1, rejected1)
plt.title(f"proposal sigma: {sigma1}")

plt.subplot(312)
plot_trace(posterior2, rejected2)
plt.title(f"proposal sigma: {sigma2}")

plt.subplot(313)
plot_trace(posterior3, rejected3)
plt.title(f"proposal sigma: {sigma3}")

plt.tight_layout()
plt.savefig("randomwalk_MH_trace.png")

Abnormally Distributed

統計解析担当のライフサイエンス研究者 -> データサイエンティスト@コンサル