tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

「Probabilistic Programming and Bayesian Methods for Hackers」を読んだ - (3章)

第3章はMCMCの原理について。

MCMCとは

詳しくはマルコフ連鎖モンテカルロ法 - Wikipediaを見るか、ベイズ統計と統計物理を参照。

MCMCマルコフ連鎖モンテカルロ法(Markov chain Monte Carlo methods)の略称で、求めたい確率分布に従って動くマルコフ過程によってその分布をサンプリングするアルゴリズム。色々細かく違いがあるアルゴリズムがあるがそれらを総称している(厳密にはマルコフ過程でないものも含む)。具体的なアルゴリズムは簡単で、

  • どこか適当にスタート地点を決める。
  • 次の移動地点を決める。
  • 分布から計算した確率でその点に移動するか决める。
  • 繰り返す

だけである。

MCMC以外の分布近似アルゴリズム

があるらしい

例: 混合分布の教師なし分類

モデリングとパラメータの推定

f:id:kuyata:20140820153159p:plain

上のような分布があるとき、この分布は2つのピークがあるように見える。これが2つの正規分布を混合したものと考えてそれぞれを推定することを考える。この分布が

  • 確率pで分布1から1-pで分布2からサンプリングする。
  • それぞれの分布は平均 \mu_{i}標準偏差 σ_{i}とする。
  • pの事前分布は[0, 1]の一様分布。
  •  \mu_{i}はそれぞれ120、190あたりがもっともらしいのでそこを中心とする分散の逆数が0.01の正規分布( τ = \dfrac{1}{σ^{2}})。
  •  τ σが[0,100]の一様分布とする。

と考えると、

p = pm.Uniform("p", 0, 1)
# 0か1を指定した確率に従ってランダムにサンプリングする
assignment = pm.Categorical("assignment", [p, 1 - p], size=data.shape[0])
centers = pm.Normal("centers", [120, 190], [0.01, 0.01], size=2)
taus = 1 / pm.Uniform("stds", 0, 1, size=2) ** 2

@pm.deterministic
def center_i(assignment=assignment, centers=centers):
    return centers[assignment]
@pm.deterministic
def tau_i(assignment=assignment, taus=taus):
    return taus[assignment]

observation = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

model = pm.Model([p, assignment, taus, centers])
mcmc = pm.MCMC(model)
mcmc.sample(50000)
# 結果を見て追加のサンプリング
mcmc.sample(100000)

# 2回めのサンプル
centers_trace = mcmc.trace("centers", chain=1)[:]
# 1回めのサンプル
prev_center_trace = mcmc.trace("centers", chain=0)[:]

MCMC.traceはchainを引数に取る。chainはデフォルトが-1でこのとき最後にMCMC.sampleを読んだ時の結果を返す。  \mu σ最尤推定はtraceしたものの平均を取ればいい。

f:id:kuyata:20140820164054p:plain

収束判定

  • 自己相関:  R(k) = Corr(x_{t}, x_{t + k})
  • MCMCは原理的にサンプリングするときに自己相関があるサンプルを取ってくる。この自己相関が高すぎる場合に分布が行ったり来たりを繰り返す状態になって収束しないことがある。
  • nサンプルに1つしかサンプルとして採用しないことで自己相関を抑えることができる。この場合より多くのサンプリングが必要になりトレードオフが生じる。
  • pymc.Matplot.plotを使うことで自己相関をプロットして視覚的に見ることができる。
from pymc.Matplot import plot as mcplot
mcplot(mcmc.trace("param"), commmnon_scale=False)

f:id:kuyata:20140820173936p:plain

左上がtrace、左下が自己相関、右が分布をプロットしたもの。

MAP(maximum a prior)を使う

MCMC.sampleだけを利用した場合、初期値の位置がうまく選択できていない場合burn in期間を長く取ることになる。初期値を最尤値に近いところで自動的にうまくとってくれるのがMAPである。

# pm.MAP(model).fit()でも良さそう
map_ = pm.MAP(model)
map_.fit()

補足

  • fitにはアルゴリズムが複数あるらしい。デフォルトはscipyのfminで、その他にはPowellの方法fmin_powellがあるらしくfit(method='fmin_powell')として指定できる。
  • 自分でvalueの値を指定してやってもいい。
  • 事前分布をうまく選ぶことが必要。収束しない場合分布の選び方が悪いことが多い。