「Pythonによるベイズ統計モデリング」まとめ #2

パラメータ $x_i$ の初期値を決める
新しいパラメータの値 $x_{i+1}$ を決め，サンプリングが簡単な分布（例：正規分布） $Q(x_{i+1}|x_i)$ からサンプルを得る。
メトロポリス-ヘイスティングスの基準

$p_a(x_{i+1}|x_i) =min \Bigl( 1, \frac{p(x_{i+1})q(x_i|x_{i+1})}{p(x_i)q(x_{i+1}|x_i)}\Bigl)$

を使って，新しいパラメータの値を受け入れる確率を計算する。
3.で計算された確率が区間[0,1]の一様分布から得られる値より大きい場合，新しい状態を受容して，そうでなければ古い状態に留まる。
1.から4.を繰り返す。

最終的に，サンプルチェーン（サンプルトレース）という数値リストを得る。これは事後分布の近似である。

Pymc3を使う

インストールについてはドキュメントを参照。

PyMC3 Documentation — PyMC3 3.6 documentation

#1の例を用いて手順を説明する。

(1)　5回投げて2回表が出た，事前分布は $Beta(1,1)$ とする

(2)　さらに5回投げて1回表が出た，事前分布は(1)の確率分布とする

(1)

まずデータを生成する。dataに「5回投げて2回表が出た」ことを記述させている。

import pymc3 as pm
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')

#(1)
data = np.array([0,1,0,0,1]) #5回投げて2回表

次にモデルを記述する。ここでは

$\theta \sim Beta(1,1)$ ， $y \sim B(n=1,p=\theta)$

であることを示す。

with pm.Model() as our_first_model:
    theta = pm.Beta('theta', alpha=1, beta=1) #事前分布，ベータ分布
    y = pm.Bernoulli('y', p=theta, observed=data) #事後分布，ベルヌーイ分布

推論ボタンを押す。

    start = pm.find_MAP() #初期値を与える（任意）
    step = pm.Metropolis() #MCMC法を用いる（任意）
    trace = pm.sample(5000, step=step, start=start) #ステップ数

これでサンプリングが始められる。

その結果を表示する。

chain = trace[burnin:] 
pm.traceplot(chain, lines={'theta':0.5}); #グラフ作成

plt.figure()

f:id:ryosuke_okubo:20190301135606p:plain

左はカーネル密度推定のグラフ（事後分布の近似），右は個々のサンプル値を表している。右のグラフは規則性がないのが望ましい。

要約量はsummaryで得られる。

pm.summary(chain)

f:id:ryosuke_okubo:20190301135845p:plain

最後に事後分布を要約する。

pm.plot_posterior(chain, kde_plot=True, ref_val=0.5)

plt.figure()

f:id:ryosuke_okubo:20190301140032p:plain
ref_val=0.5は，参照点が0.5であることを示す。

(2)も同じようにやってみるとよい。

第3章　複数パラメータの取り扱いと階層モデル

基本的には第2章の延長線として読めるので，ここでは概要のみとする。

正規性の仮定

先ほどの例ではベータ分布と2項分布によるモデルを用いた。この他，正規分布からなるモデルが構築できる。

$\mu \sim Uniform(l,h)$ ， $\sigma \sim HalfNormal(\sigma_{\sigma})$

$y \sim Normal(\mu,\sigma)$

もし外れ値が気になるのであれば，Studentのt分布を用いる。

$\mu \sim Uniform(l,h)$ ， $\sigma \sim HalfNormal(\sigma_{\sigma})$

$\nu \sim Exponential(\lambda)$

$y \sim StudentT(\mu,\sigma.\nu)$

Studentのt分布では頑健推定ができる。

グループ間比較

ベイジアン分析においてはp値は用いず，代わりに効果量を用いる。効果量を測定する方法として，コーエンのdや優越率などが用いられる。

階層モデル

モデルのパラメータを入れ子にすることで，階層モデルが構築できる。本書の例は以下の通りである。

$\alpha \sim HalfCauchy(\beta_{\alpha})$ ， $\beta \sim HalfCauchy(\beta_{\beta})$

$\theta \sim Beta(\alpha,\beta)$

$y \sim Bern(\theta)$

参考：はてなブログにコードを挿入する方法

【はてなブログ】ソースコードを貼り付ける方法！（言語名＆行番号表示） - 黒木ノ水岩

次回↓

ryosuke-okubo.hatenablog.com

十の並列した脳

何でも勉強する，毎週月木曜に投稿予定

「Pythonによるベイズ統計モデリング」まとめ #2

第2章　確率プログラミング－PyMC3入門

推論エンジン

Pymc3を使う

第3章　複数パラメータの取り扱いと階層モデル

正規性の仮定

グループ間比較

階層モデル

第2章 確率プログラミング－PyMC3入門

推論エンジン

Pymc3を使う

第3章 複数パラメータの取り扱いと階層モデル

正規性の仮定

グループ間比較

階層モデル

第2章　確率プログラミング－PyMC3入門

第3章　複数パラメータの取り扱いと階層モデル