So-net無料ブログ作成

PythonでMCMC(メトロポリス法)-ベイズ統計学 [統計学]

最近は状態空間モデリングの勉強をしています。その一環できちんとMCMCを理解しようとあれこれ資料を漁っていました。マルコフ連鎖モンテカルロ法を実装してみようというブログ記事や「計算統計 2 マルコフ連鎖モンテカルロ法とその周辺」などを読みつつ、メトロポリス法というシンプルなMCMCをPythonで実装してみました。

■計算条件
・2変量正規分布からサンプリング
【2変量正規分布】


・b = 0.5 (ここでは適当にこの値にした)
・burn-in = 0
-サンプリングの挙動を確認するためバーンインは0とした

■Pythonコード
#!/usr/bin/env python
#-*- coding:utf-8 -*-

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

def metropolis(func, start=10.0, delta=2.0, burn_in=0, seed=1):
    metro_iter = _metropolis(func, start, delta, seed)
    for i, value in enumerate(metro_iter):
        if i > burn_in:
            yield value
            break
    for value in metro_iter:
        yield value
       
def _metropolis(func, start, delta, seed):
    random.seed(seed)
    initial_x1 = start
    initial_x2 = start
    initial_fx = func(initial_x1, initial_x2)
    proposed_x1, proposed_x2, proposed_fx = 0.0, 0.0, 0.0
    for j in range(1,1000):
        proposed_x1 = initial_x1 + random.uniform(-delta, delta)
        proposed_x2 = initial_x2 + random.uniform(-delta, delta)
        proposed_fx = func(proposed_x1, proposed_x2)
        ratio = proposed_fx / initial_fx
        if ratio > 1.0 or ratio > random.random():
            initial_x1, initial_x2, initial_fx = proposed_x1, proposed_x2, proposed_fx
            yield proposed_x1, proposed_x2, proposed_fx
        else:
            yield initial_x1, initial_x2, initial_fx

def main():
    b = 0.5
    Z = 2*math.pi/math.sqrt(1.0-b**2.0)
    f = lambda x1, x2: math.exp(-(x1**2-2*b*x1*x2+x2**2)/2)/Z
    xb1 = []
    xb2 = []
    fxb1 = []
    for x1, x2, fx in metropolis(f):
        xb1.append(x1)
        xb2.append(x2)
        fxb1.append(fx)
    
    tx = range(1,len(fxb1)+1)
    md = [np.median(fxb1)]*(len(fxb1)) 
    
    plt.plot(xb1,xb2,"bo")
    plt.axis([-10, 10, -10, 10])
    plt.show()
    plt.plot(tx, fxb1)
    plt.plot(tx, md, "r")
    plt.axis([0, 1000, 0, 0.14])  
    plt.show()

if __name__ == '__main__':
    main()


■結果
・2変量の動き
figure_1.png
#初期値から数サンプルはふらふらと彷徨うが、しばらくすると最大値付近に収束

・p(x1,x2)の収束状況
figure_2.png
#時系列にグラフ化すると収束性を確認することができる

■感想
・メトロポリス法はとシンプルなアルゴリズムで理解しやすい
・この調子でメトロポリスヘイスティング法とギブスサンプラー法も理解する予定
・GAなんかを併用する手法はないのかな?

■参考文献

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)

  • 作者: 伊庭 幸人
  • 出版社/メーカー: 岩波書店
  • 発売日: 2005/10/27
  • メディア: 単行本


・伊庭先生の解説が一番分かりやすいと思う。これが少し難易度が高いなら久保先生のベイズ本などを最初に読むといいです

ベイズ統計データ解析 (Rで学ぶデータサイエンス 3)

ベイズ統計データ解析 (Rで学ぶデータサイエンス 3)

  • 作者: 姜 興起
  • 出版社/メーカー: 共立出版
  • 発売日: 2010/07/24
  • メディア: 単行本


・あんまり注目されていないような気がするんだけど、この本も名著です。MCMCのRコードもあります。状態空間モデルのRコードはとても重宝しています。北川先生の時系列本と合わせて読むと理解がより深まりますね!



タグ:MCMC ベイズ
nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

トラックバックの受付は締め切りました
メッセージを送る

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。

×

この広告は1年以上新しい記事の更新がないブログに表示されております。