Pythonで学ぶ！ ベイズ統計 第6回

• 前回 : Pythonで学ぶ！ ベイズ統計 第5回 - ma-hi-ro Hobby

Pythonで学ぶ！ ベイズ統計 第6回

前回は、実計算で確率分布のベイズ更新を行い、パラメータ推定を行った
今回は、MCMCを用いて、同様にパラメータの推定を行う

PyMC3のインストール

普通にpipでインストールできた

pip install pymc3

おまけ Python環境について

個人的に、Python環境は下記を推奨

• Python公式からインストール
• venvで仮想環境構築
• 基本pipのみでモジュール管理

交通事故の例1

前回の交差点での交通事故の例をMCMCで求める

まずは、PyMC3モジュールのインポート

import pymc3 as pm
import warnings
warnings.simplefilter('ignore', FutureWarning)

モデルの設定と推定
今回は、1データ（4）のみ

data = [4]
with pm.Model() as model:
    lambda_ = pm.Gamma('lambda', alpha=1.6, beta=1)
    pm.Poisson('x', mu=lambda_, observed=data)
    trace = pm.sample(random_seed=1)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 32 seconds.

デフォルトの設定だと、調整に1000回、推定に1000回乱数を生成して、それを４サイクル実行している

表示

pm.plot_trace(trace, figsize=(12, 5))

4サイクル分のの推定値
左が分布で、右が1000回の変遷

結果のSummaryの表示

pm.summary(trace)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
lambda	2.771	1.161	0.774	4.789	0.027	0.019	1735.0	2465.0	1.0

前回の計算結果だと、平均 2.8 / 分散 1.4 に対して、2.771 / 1.348 (=1.161x1.161)

交通事故の例2

データ数を増やした場合

data = [4, 3, 1, 2, 1, 5]
with pm.Model() as model:
    lambda_ = pm.Gamma('lambda', alpha=1.6, beta=1)
    pm.Poisson('x', mu=lambda_, observed=data)
    trace = pm.sample(random_seed=1)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 31 seconds.

pm.plot_trace(trace, figsize=(12, 5))

結果のSummaryの表示

pm.summary(trace)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
lambda	2.509	0.602	1.407	3.613	0.014	0.01	1715.0	2168.0	1.0

前回の計算結果だと、平均 2.51 / 分散 0.36 に対して、2.509 / 0.362 (=0.602x0.602)

参考文献

この記事は以下の書籍を参考にしましたが、
私の拙い知識で書いておりますので、誤り等ありましたらご指摘ください

Pythonで動かして学ぶ！あたらしいベイズ統計の教科書