Gibbs Sampling - whats-mind @ ウィキ

最近Dirichlet Processというのを良く聞きます。
クラスタリングというとあらかじめのクラスタ数Kで決める手法が多いのですが、
この手法はKそのものを推定する手法だそうです。すごいですね～
ということで、Dirichlet Processの勉強をしたくて
その学習手法の中で出てくるギブスサンプリングを
簡単な例でpythonで実装してみました。

ギブスサンプリングを用いてX～N(μ,σ^2)から
出力されたサンプルからμとσを推定します。
μの事前分布をガウス分布、σの事前分布を逆ガンマ分布としています。

#!/usr/bin/python
#coding:utf-8
from numpy import *
from scipy import ndimage
from scipy.stats import norm
from scipy.stats import invgamma

# サンプル数
n_sample = 30

# ステップ数
n_step = 1000

# 一次元ガウス分布の乱数
def r_gaussian(mu, sigma2, size=0):
    r_norm = norm.rvs(size=size)
    return [mu + sqrt(sigma2)*r for r in r_norm]

# 逆ガンマ分布
def r_invgamma(alpha, beta, size=0):
    r_invgamma0 = invgamma.rvs(alpha, size=size)
    return [beta*r for r in r_invgamma0]

# 推定したい平均と分散の2乗
mu_c   = 3.0
sig2_c = 1.0
sample = r_gaussian(mu_c, sig2_c, n_sample)

# muの事前分布パラメタ
mu_o   = 5.0
tau2_o = 10.0

# sigma2の事前分布のパラメタ
nu_o  = 10.0
lmd_o = 30.0

# 推定値
mu_e   = 0.0
sig2_e = 1.0

avg_sample  = mean(sample)
sig2_sample = ndimage.variance(sample)

data_mu   = []
data_sig2 = []
# Gibbs Sampling
for n in range(n_step):
    mu_e = r_gaussian((n_sample*avg_sample/sig2_e+mu_o/tau2_o)/(n_sample/sig2_e+1.0/tau2_o), 1.0/(n_sample/sig2_e+1.0/tau2_o), 1)[0]
    sig2_e = r_invgamma((n_sample+nu_o)/2.0, (sig2_sample+n_sample*(mu_e-avg_sample)**2+lmd_o)/2.0, 1)[0]
    data_mu.append(mu_e)
    data_sig2.append(sig2_e)
#描画
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
#ax.plot(sample)
ax.plot(data_mu)
ax.plot(data_sig2)
ax.set_ylim([0.0, 5.0])
plt.show()

緑が分散で青が平均の推定値です。