【M-H法】Python で MCMC を書いてみた【Gibbs-Sampler】

| 1件のコメント

この記事は Python Advent Calendar 2015 19日目の記事です。
MC法 (Monte Carlo Method) からスタートして, MCMC (Markov Chain Monte Carlo)を Python で書いてみます。

MC法による円周率の近似計算

MC法 (Monte Carlo Method)は乱数を利用した計算アルゴリズム。
MC法の hello world 的な例題として, 円周率の近似値を計算する。

import numpy as np
import matplotlib.pyplot as plt

N = 100000

x = np.random.uniform(-1.0, 1.0, N)
y = np.random.uniform(-1.0, 1.0, N)

l = []

for i in xrange(N):
    if (x[i]**2 + y[i]**2) < 1.0:
        l.append(True)
    else:
        l.append(False)

print 2.0**2 * float(l.count(True)) / float(N) # 3.13632

plt.scatter(x, y, c=l, s=5, edgecolor='None')
plt.xlim(-1.0, 1.0)
plt.ylim(-1.0, 1.0)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

正方形の領域に一様分布から 10万点を取り出し打ち出す。
円の中に入った点の個数をカウントすると正方形の面積と円の面積の比が求まるので, これに正方形の面積を掛けると, 円の面積が出てくる。

monte-carlo-2dim

結果は 3.13632 と求まり, 半径が 1 なので円周率(=面積)の近似値となる。
同様に球の体積でも MC法で 4.19312 と近似値計算は上手くいっている。

monte-carlo-3dim

MC法による超球の体積近似計算

円や球ではMC法でそれなりの精度で近似値が得られた。
同様の方法で, 次元を上げていった時の n次元超球の体積が上手く求まるか試してみる。

半径 r の n次元超球の体積は下記で求まる。

     \begin{eqnarray*}   V_{n}(r) = \frac{\pi^\frac{n}{2}}{\Gamma\frac{n}{2}+1}r^n \end{eqnarray*}

超球の体積をMC法で 1次元から24次元まで求め, 理論値と比較してみる。

import math
import numpy as np
import matplotlib.pyplot as plt

def V(n):
    return math.pi**(n/2.0) / math.gamma(n/2.0+1.0)

def count_point(n):
    x = []
    count = 0

    for i in xrange(n):
        x.append(np.random.uniform(-1.0, 1.0, N))

    for i in xrange(N):
        r = 0.0
        for j in xrange(n):
            r += x[j][i]**2.0
        if r < 1.0:
            count += 1

    return count

tv = []
mc = []
N = 100000

for n in xrange(1,25):
    c = count_point(n)
    tv.append(V(n))
    mc.append(2.0**n * float(c) / float(N))

    print n, 'Dim'
    print 'Theoretical Value:', V(n)
    print 'Monte Carlo :', 2.0**n * float(c) / float(N)

x = np.arange(1, 25, 1)
plt.plot(x, tv)
plt.plot(x, mc)
plt.xlabel('dim')
plt.ylabel('V(r=1)')
plt.grid(True)
plt.show()

下記が結果で, MC法では 15次元で体積が求まらなくなってしまった。

monte-carlo-high-dim

このように, 高次元になるほど単純な MC法では超球の内部へ稀にしか入らず, 近似計算が上手くいかなくなる。(次元の呪い)

MCMC

MCMC (Markov Chain Monte Carlo, マルコフ連鎖モンテカルロ法) はMC法の中でも, 調べたい分布を不変分布 [3]とするマルコフ連鎖 [4]をシミュレートし, サンプルを生成する方法である。生成したサンプルの分布を調べることで, パラメータ推定に使うこともできる。

確率密度の大きいところを探しながらサンプルを生成できる点が特徴で, 単純な MC法によるn次元超球の体積計算のように無駄なサンプルが増えてしまう影響を抑制できる。

MCMCに関する性質・条件 (既約性, エルゴード性, 詳細釣り合い条件) 等については専門書を参考に。

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

メトロポリス・ヘイスティングス法 (Metropolis-Hastings algorithm, M-H法) は ギブス・サンプラー (Gibbs sampling, Gibbs sampler)と メトロポリス法 (Metropolis method)を含む MCMCの代表的な手法のひとつ。

ギブス・サンプラー

ギブス・サンプラーは新しいサンプルを全て受容できる効率的なサンプリング手法。

例として二変量正規分布で試してみる。(数式引用 [2])

     \begin{eqnarray*}   p(x_1, x_2) = \frac{1}{Z} exp(-\frac{x_1^2 - 2bx_1x_2 + x_2^2}{2}) \end{eqnarray*}

二変量正規分布を3Dプロットしてみる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

def range_ex(start, end, step):
    while start + step < end:
        yield start
        start += step

# P(x) : Target distribution
def P(x1, x2, b):
    return np.exp(-1/2 * (x1**2 - 2*b*x1*x2 + x2**2))

xs = []
ys = []
zs = []
b = 0.5

for i in range_ex(-3, 3, 0.1):
    for j in range_ex(-3, 3, 0.1):
        xs.append(i)
        ys.append(j)
        zs.append(P(i, j, b))

ax = Axes3D(plt.figure())
ax.scatter3D(xs, ys, zs, s=3, edgecolor='None')
plt.show()

bivariate-normal-distribution

上記の二変量正規分布の x_1, x_2について, 一方を固定した場合の条件付き確率がわかっているとすると, ギブス・サンプラーは小さい N の大きさで効率的にサンプルを得られる。一方を固定しているので, 平面上をカクカク動くのが特徴。

def gibbs(N, thin):
    s = []
    x1 = 0.0
    x2 = 0.0
    for i in range(N):
        for j in range(thin):
            x1 = np.random.normal(b * x2, 1) # P(x1|x2)
            x2 = np.random.normal(b * x1, 1) # P(x2|x1)
        s.append((x1,x2))
    return np.array(s)

N = 3000
thin = 500
burn_in = 0.2

sample = gibbs(N, thin)

plt.scatter(
    sample[int(N * burn_in):,0],
    sample[int(N * burn_in):,1],
    s=3,
    edgecolor='None'
)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()

print 'x1 :', np.mean(sample[:,0]), np.var(sample[:,0]) # x1 : -0.0176636194179 1.32838943108
print 'x2 :', np.mean(sample[:,1]), np.var(sample[:,1]) # x2 : -0.0190149176882 1.28633687028

初期状態の影響が残っている最初の方の列を捨てると, 残りの列は定常分布, つまり解析したい分布からのサンプル列と考えられる。今回は, burn_in = 0.2 として最初の 20%は捨てている。

gibbs-sampling-normal-distribution

メトロポリス法

ギブス・サンプラーでは条件付確率を用いたが, 条件付確率を求めるのが容易でない場合は, メトロポリス法を用いる。

条件付確率の代わりに目標分布 p(z)を近似する簡単な提案分布 q(z) を用いる。ただし, q(z)から得られた 新しいサンプルの候補について, 受容か棄却かを決定するのは p(z)を使って求める。

遷移する確率は a が 1 より大きいと必ず遷移し, 1以下の場合はその確率 ([0,1] の一様分布からの乱数と比較した結果)で遷移する。

import copy
import numpy as np
import matplotlib.pyplot as plt

# P(x) : Target distribution
def P(x1, x2, b):
    return np.exp(-0.5 * (x1**2 - 2*b*x1*x2 + x2**2))

# Q(x) : Proposal distribution
def Q(c, mu1, mu2, sigma):
    return (c[0] + np.random.normal(mu1, sigma), c[1] + np.random.normal(mu2, sigma))

def metropolis(N):
    current = (10, 10)
    sample = []
    sample.append(current)
    accept_ratio = []

    for i in range(N):
        candidate = Q(current, mu1, mu2, sigma)

        T_prev = P(current[0], current[1], b)
        T_next = P(candidate[0], candidate[1], b)
        a = T_next / T_prev

        if a > 1 or a > np.random.uniform(0, 1):
            # Update state
            current = copy.copy(candidate)
            sample.append(current)
            accept_ratio.append(i)

    print 'Accept ratio :', float(len(accept_ratio)) / N
    return np.array(sample)

b = 0.5
mu1 = 0
mu2 = 0
sigma = 1

N = 30000
burn_in = 0.2

sample = metropolis(N)

plt.scatter(
    sample[int(len(sample) * burn_in):, 0],
    sample[int(len(sample) * burn_in):, 1],
    alpha=0.3,
    s=5,
    edgecolor='None'
)
plt.title('MCMC (Metropolis)')
plt.show()

fig = plt.figure(figsize=(15, 6))

ax = fig.add_subplot(121)
plt.hist(sample[int(N * burn_in):,0], bins=30)
plt.title('x')

ax = fig.add_subplot(122)
plt.hist(sample[int(N * burn_in):,1], bins=30)
plt.title('y')
plt.show()

print 'x :', np.mean(sample[:,0]), np.var(sample[:,0]) # x : 0.02767180711 1.36504061679
print 'y :', np.mean(sample[:,1]), np.var(sample[:,1]) # y : 0.0309528474353 1.37879772867

新しいサンプルの候補が, 受容された比率は, 0.5623 となった。

metropolis-algorithm-normal-distribution

MCMCで生成した多数のサンプルの分布から, 平均・分散といった要約統計量やEAP推定量, ベイズ信頼区間といった欲しい情報を求めることができる。

mcmc-result-distribution

Codeは Githubに置いた。

おわりに

間違いを見つけられましたら, ご指摘いただけると嬉しいです。最後に, 参考にさせて頂いた書籍を何冊か載せておきます。

MCMC以外にも, オンライン学習や正則化の辺りが分かり易かったです。

MCMC (動的モンテカルロ法)は元々, 統計物理との関連が深いことをこの書籍で知りました。

ベイズ統計の歴史的背景や, Stanでも実装されている ハミルトニアン・モンテカルロ (HMC)の話が分かり易かったです。


[1] Gibbs sampler in various languages (revisited)
[2] 第11章 サンプリング法
[3] 不変分布とは, (ベイズ統計的な意味で) ある一定の確率分布へと収束すること。
[4] t の状態が t-1 からのみの条件付確率で決定することをマルコフ性と言い, その確率過程をマルコフ過程と言う。その中でも離散的な状態または時間の場合がマルコフ連鎖。
[5] 例えば, 正規分布の場合, MCMCで得られた統計量から確率密度を確認したければ R で curve(dnorm(x, mu, sigma))) とかでできる (あまり意味はないけど...)
[6] マルコフ連鎖モンテカルロ法の最近の展開 (2001)

1件のコメント

  1. ピンバック: 2015年を振り返る | FiS Project

コメントを残す

Required fields are marked *.