StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第11章「練習問題」
第11章「離散値をとるパラメータを使う」
書籍の著者 松浦健太郎 先生
この記事は、テキスト第11章「離散値をとるパラメータを使う」の「練習問題」の PyMC5写経 を取り扱います。
なお、次の練習問題は写経を省略いたしました。
練習問題(1)：Stan特有の文法
練習問題(4)：数式による証明
練習問題(5)：Stan特有の文法
練習問題(6)：PyMCの工夫点を思いつかない
はじめに
StanとRでベイズ統計モデリングの紹介
この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。
テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！
テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
PyMC環境の準備
Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。
11章 練習問題
インポート
### インポート
# 数値・確率計算
import pandas as pd
import numpy as np
# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')
練習問題(2)
11.1.2項のベルヌーイ分布の例題のデータを読み込みます。
### データの読み込み ◆データファイル11.1 data-coin.txt
# Y: 回答（1:Yes） サイコロを投げて1：喫煙経験を答える、その他の目：常にYesを答える
data1 = pd.read_csv('./data/data-coin.txt')
print('data1.shape: ', data1.shape)
display(data1.head())【実行結果】
モデルの定義です。
### モデルの定義 （周辺化消去しない）
with pm.Model() as model1:
### データ関連定義
## coordの定義
model1.add_coord('data', values=data1.index, mutable=True)
## dataの定義
# 目的変数 Y
Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
### 事前分布
# サイコロ投げの結果 1以外の目:0, 1の目:1
dice1 = pm.Bernoulli('dice1', p=1/6, dims='data')
# 喫煙確率 q
q = pm.Uniform('q', lower=0, upper=1)
# Yesと回答する確率 θ サイコロの目が1以外:1, 1の目:q
theta = pt.stack([1, q])
### 尤度関数
obs = pm.Bernoulli('obs', p=theta[dice1], observed=Y, dims='data')モデルの定義内容を見ます。
### モデルの表示
model1【実行結果】
### モデルの可視化
pm.model_to_graphviz(model1)【実行結果】
MCMCを実行します。
### 事後分布からのサンプリング 1分
# dice1の初期値設定
initvals = {'dice1': np.ones(len(data1))}
# MCMCの実行 ※PyMC標準のNUTSサンプラーを使用
with model1:
idata1 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
initvals=initvals, random_seed=1234)【実行結果】
事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。
### r_hat>1.1の確認
# 設定
idata_in = idata1 # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
収束条件を満たしています。
事後統計量を表示します。
### 推論データの要約統計情報の表示
pm.summary(idata1, hdi_prob=0.95, round_to=3)【実行結果】
トレースプロットを描画します。
### トレースプロットの表示
pm.plot_trace(idata1, compact=True)
plt.tight_layout();【実行結果】
《解答》
ｑの事後統計量を算出します。
事後統計量算出関数の定義です。
### median, 2.5%, 97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
probs = [50, 2.5, 97.5]
columns = ['median', '2.5%', '97.5%']
quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
quantiles.columns=columns
return quantiles事後統計量を算出します。
### qの中央値・95%信用区間の算出 ◆解答
vars = ['q']
param_samples = idata1.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(3))【実行結果】
練習問題(3) A)
２枚のコインのうち１枚は表が出る確率を 0.5 とし、もう１枚は表が出る確率を 0.4 とするケースです。
11.1.2項のベルヌーイ分布の例題のデータを用います。
モデルの定義です。
### モデルの定義 （周辺化消去しない）
with pm.Model() as model2:
### データ関連定義
## coordの定義
model2.add_coord('data', values=data1.index, mutable=True)
## dataの定義
# 目的変数 Y
Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
### 事前分布
# サイコロ投げの結果 ２回coinトスして、2回表:1 , それ以外:0
# ２枚表が出る確率＝0.5 * 0.4 = 0.2
coin2 = pm.Bernoulli('coin2', p=0.2, dims='data')
# 喫煙確率 q
q = pm.Uniform('q', lower=0, upper=1)
# Yesと回答する確率 θ 右以外:1, 2回表:q
theta = pt.stack([1, q])
### 尤度関数
obs = pm.Bernoulli('obs', p=theta[coin2], observed=Y, dims='data')モデルの定義内容を見ます。
### モデルの表示
model2【実行結果】
### モデルの可視化
pm.model_to_graphviz(model2)【実行結果】
MCMCを実行します。
### 事後分布からのサンプリング 1分
# coin2の初期値設定
initvals = {'coin2': np.ones(len(data1))}
# MCMCの実行 ※PyMC標準のNUTSサンプラーを使用
with model2:
idata2 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
initvals=initvals,
random_seed=1234)【実行結果】
事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。
### r_hat>1.1の確認
# 設定
idata_in = idata2 # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
収束条件を満たしています。
事後統計量を表示します。
### 推論データの要約統計情報の表示
pm.summary(idata2, hdi_prob=0.95, round_to=3)【実行結果】
トレースプロットを描画します。
### トレースプロットの表示
pm.plot_trace(idata2, compact=True)
plt.tight_layout();【実行結果】
《解答》
ｑの事後統計量を算出します。
### qの中央値・95%信用区間の算出 ◆解答
vars = ['q']
param_samples = idata2.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(3))【実行結果】
練習問題(3) B)
２枚のコインのうち１枚は表が出る確率を 0.5 とし、もう１枚は表が出る確率を p_coin（無情報事前分布） とするケースです。
11.1.2項のベルヌーイ分布の例題のデータを用います。
モデルの定義です。
### モデルの定義 （周辺化消去しない）
with pm.Model() as model3:
### データ関連定義
## coordの定義
model3.add_coord('data', values=data1.index, mutable=True)
## dataの定義
# 目的変数 Y
Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')
### 事前分布
# サイコロ投げの結果 ２回coinトスして、2回表:1 , それ以外:0
# ２枚表が出る確率＝0.5 * pCoin
pCoin = pm.Uniform('pCoin', lower=0, upper=1)
coin2 = pm.Bernoulli('coin2', p=0.5 * pCoin, dims='data')
# 喫煙確率 q
q = pm.Uniform('q', lower=0, upper=1)
# Yesと回答する確率 θ 右以外:1, 2回表:q
theta = pt.stack([1, q])
### 尤度関数
obs = pm.Bernoulli('obs', p=theta[coin2], observed=Y, dims='data')モデルの定義内容を見ます。
### モデルの表示
model3【実行結果】
### モデルの可視化
pm.model_to_graphviz(model3)【実行結果】
MCMCを実行します。
### 事後分布からのサンプリング 1分10秒
# coin2の初期値設定
initvals = {'coin2': np.ones(len(data1))}
# MCMCの実行 ※PyMC標準のNUTSサンプラーを使用
with model3:
idata3 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
initvals=initvals,
random_seed=1234)【実行結果】
事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。
### r_hat>1.1の確認
# 設定
idata_in = idata3 # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
収束条件を満たしています。
事後統計量を表示します。
### 推論データの要約統計情報の表示
pm.summary(idata3, hdi_prob=0.95, round_to=3)【実行結果】
トレースプロットを描画します。
### トレースプロットの表示
pm.plot_trace(idata3, compact=True)
plt.tight_layout();【実行結果】
《解答》
ｑの事後統計量を算出します。
### qの中央値・95%信用区間の算出 ◆解答
vars = ['q']
param_samples = idata3.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(3))【実行結果】
第11章 練習問題は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。