PyMC3 の使い方 (ベイズ推論)

久しぶりに使おうとすると忘れるのでメモしておきます。

PyMC3 とは？

• Python で使える確率的プログラミングのライブラリ
• ベイズ推論に利用できる

インストール

pip install pymc3

基本的な書き方

import numpy as np
from matplotlib import pyplot as plt

# 試行回数と観測結果(回数)
N = 100
a = 5

with pm.Model() as model:
  # 事前分布
  #  範囲が 0 〜 1 の一様分布とする
  theta = pm.Uniform('theta', lower=0, upper=1)

  # 尤度関数
  #   ベルヌーイ試行を N 回行った時の観測回数 a が従う分布として二項分布を設定する
  obs = pm.Binomial('a', p=theta, n=N, observed=a)

  # 推論を実行し、事後分布から 5000*2 サンプルを得る
  trace = pm.sample(5000, chains=2)

• pm.sample では MCMC でサンプルを得て trace に結果が格納されている

サンプル(事後分布)の可視化

with model:
  pm.traceplot(trace)

• 左側の図には対象の確率変数について推定された事後分布が表示される
• 右側にはバーンイン(burn-in)以降に得られたサンプルの軌跡が表示される
  • バーンイン(burn-in): 探索の初期値の影響を取り除くためにサンプルを捨ててしまう期間のこと

要約統計量の表示

with model:
  # HDI(highest density interval) として 95% HDI を設定
  print(pm.summary(trace, hdi_prob=0.95)

要約統計量を付きでサンプル(事後分布)を可視化

with model:
  pm.plot_posterior(trace, hdi_prob=0.95)

参考資料

• PyMC3 Documentation