ギブスサンプリングの実装

長いことわかっていなかったのだが、実装してようやく頭に入って来た気がする。

今回は「道具としてのベイズ統計」涌井良幸　日本実業社
の説明にあるギブスサンプリングをPythonで実装してみた。
パラメータは全部、本にある通り。

#! /usr/bin/python                                                                                                                                                           
# -*- coding:utf-8 -*-                                                                                                                                                       
"""                                                                                                                                                                          
Gibbs samplingの実装。「道具としてのベイズ統計」のパラメータをそのまま流用。                                                                                     
"""
import sys, math, random
import numpy as np
import pylab as pl
import scipy.stats as stats

number_of_data=30
Q=138.22
average=5.11

n_0=0.02
S_0=1.0
mu_0=5.0
m_0=0.25

mu_1=(number_of_data*average+(m_0*mu_0))/(m_0+number_of_data)
m_1=30.25
n_1=float(number_of_data)
n_1_times_S_1=138.24

def sample_mu_from_normal_distri_of_after(sigma):
    mu_in_mu_distribution_after=mu_1
    sigma_in_mu_distribution_after=(sigma)/m_1
    mu_now=stats.norm.rvs(mu_in_mu_distribution_after, sigma_in_mu_distribution_after)
    return mu_now
def sample_sigma_from_inverse_gamma_distri_of_after(mu_now):
    mu_in_sigma_distribution_after=(n_1+1.0)/2.0
    sigma_in_sigma_distribution_after=((n_1_times_S_1)+m_1*((mu_now-mu_1)**2))/2.0
    #ちょっとこの関数には妙な感じがするので、別の方法で                                                                                                                      
    #sigma_now=stats.invgamma.rvs(sigma_in_sigma_distribution_after, loc=mu_in_sigma_distribution_after)                                                                     
    #sigma_now=math.sqrt(sigma_now)                                                                                                                                          
    sigma_now=(1.0/random.gammavariate(mu_in_sigma_distribution_after, 1.0/sigma_in_sigma_distribution_after))
    return sigma_now
def main(ITER_NUM):
    sigma_now=4.0
    mu_list=[]
    sigma_list=[]
    for i in range(0, ITER_NUM):
        mu_now=sample_mu_from_normal_distri_of_after(sigma_now)
        sigma_now=sample_sigma_from_inverse_gamma_distri_of_after(mu_now)
        mu_list.append(mu_now)
        sigma_list.append(sigma_now)

    del mu_list[0:999]
    del sigma_list[0:999]
    mu_summation=0
    for x in mu_list:mu_summation+=x
    average_mu=float(mu_summation)/len(mu_list)
    sigma_summation=0
    for x in sigma_list:sigma_summation+=x
    average_sigma=sigma_summation/len(sigma_list)
    print average_mu, average_sigma
if __name__=='__main__':
    main(ITER_NUM=3000)

実行すると、尤度関数の
平均5.10653864685
分散4.76013929248

となり、まあ、しっかりできてますね。というのが確認できた。
ところで、気になるのは、逆ガンマ分布のサンプリング（つまり、乱数発生のさせ方）。
引用書の中ではExcelには逆ガンマ分布の乱数発生コマンドがないので、ガンマ分布の引数にラムダ（中心）の逆数を与え、乱数値の逆数をとるように。と指示されている。

けど、PythonにはどうやらScipyに逆ガンマ分布のコマンドが用意されているようで、
http://gaezipserver.appspot.com/python/sci/scipy/generated/scipy.stats.invgamma.html
しかも、乱数発生のコマンドもある。
けど、これを使うと、どうしても値が合わなくなってしまうのだ。
なので、せっかくなのだが、使わずに引用書の指示通りにした。
（該当箇所はコメントアウトしてある）

これも全部、英語が読めないのが悪いんだ....

類似して、Javaでも実装されている方もいらっしゃるので
http://d.hatena.ne.jp/yokoken00/20100128