StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第10章「10.1.5 ウサギとカメ」
第10章「収束しない場合の対処法」
書籍の著者 松浦健太郎 先生
この記事は、テキスト第10章「収束しない場合の対処法」・10.1節「パラメータの識別可能性」の 10.1.5項「ウサギとカメ」の PyMC5写経 を取り扱います。
今回はうまくモデリングできず、おかしな結果になりました。
類題のモデル式10-4（将棋）でリベンジします🦾
テキストは第10章で 収束に向けた工夫を取り扱っています。
座学や数式モデルのみの掲載項があれば、Stanコードの掲載項もあります。
PyMC化は主にStanコードが明示されているモデル式を対象にして実施します。
はじめに
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を動かすまでの準備」章をご覧ください。
10.1.5 ウサギとカメ
Stanは変数 ordered 型を用いて変数の順序の維持することができるようです。
PyMCの類似概念には、確率変数の引数 transform で用いることができる「pymc.distributions.transforms.ordered」があります。
ただし中々の気難し屋さんなので、私はうまく使いこなすことができません。
今回のPyMC写経では、順序を無視してモデリングしました。
一応収束したものの、推論が不適切な感じになっています。
インポート
### インポート
# 数値・確率計算
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
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')
データの読み込み・確認
サンプルコードのデータを読み込みます。
### データの読み込み ◆データファイル10.2 data-usagitokame.txt
# Loser:負けたプレイヤーのインデックス, Winner:勝ったプレイヤーのインデックス
# インデックスの値 1:カメ, 2:ウサギ
data = pd.read_csv('./data/data-usagitokame.txt')
print('data.shape: ', data.shape)
display(data.head())【実行結果】
勝敗を集計してみます。
### 勝敗の集計
display(data['Winner'].value_counts().sort_index().to_frame()
.rename({1:'カメ', 2:'ウサギ'}, axis=0))【実行結果】
テキストの sim-model10-3.R に相当するシミュレーションを実行します。
### シミュレーション ◆sim-model10-3.R
## 初期値設定
# 乱数生成器の設定
rng = np.random.default_rng(seed=1234)
# ゲーム数の設定
G = 30
# パフォーマンスの平均値の設定 [カメ, ウサギ]
mu_pf = np.array([0, 1.5])
## データ作成
# パフォーマンス値の作成 1列目：カメ、2列目：ウサギ
pf = rng.normal(loc=mu_pf, scale=1, size=(G, 2))
# 勝者データの作成 カメが勝利:1, ウサギが勝利2
d = (pf[:, 0] < pf[:, 1]) + 1
# データフレーム化
d = pd.DataFrame({
'カメのパフォーマンス': pf[:, 0],
'ウサギのパフォーマンス': pf[:, 1],
'Loser': np.where(d==1, 2, 1) ,
'Winner': np.where(d==1, 1, 2) , })
# データの表示
display(d['Winner'].value_counts().sort_index().to_frame()
.rename({1:'カメ', 2:'ウサギ'}, axis=0))
display(d)【実行結果】
PyMCのモデル定義
PyMCでモデル式10-2を実装します。
モデルの定義です。
観測値を設定する確率分布が無いモデルです。
カメの強さの平均値 mukame には 0 を設定します。
### モデルの定義 ◆モデル式10-3 model10-3.stan
with pm.Model() as model:
### データ関連定義
## coordの定義
# データのインデックス
model.add_coord('game', values=data.index, mutable=True)
# カメ・ウサギ
model.add_coord('kameUsagi', values=['カメ', 'ウサギ'], mutable=True)
# 負け・勝ち
model.add_coord('losewin', values=['loser', 'winner'], mutable=True)
## dataの定義
# loser 0始まり化 0:カメ, 1:ウサギ
loser = pm.ConstantData('loser', value=data['Loser'].values - 1,
dims='game')
# winner 0始まり化 0:カメ, 1:ウサギ
winner = pm.ConstantData('winner', value=data['Winner'].values - 1,
dims='game')
### mu
muKame = 0
muUsagi = pm.Uniform('muUsagi', lower=-10, upper=10)
mu = pm.Deterministic('mu', pt.stack([muKame, muUsagi]), dims='kameUsagi')
### パフォーマンス ※観測値が無い
performance = pm.Normal('performance',
mu=pt.stack([mu[loser], mu[winner]]).T,
sigma=1,
dims=('game', 'losewin'),
# orderedを含めるとMCMC実行時にAssertionErrorが発生
# transform=pm.distributions.transforms.ordered
)モデルの定義内容を見ます。
### モデルの表示
model【実行結果】
### モデルの可視化
pm.model_to_graphviz(model)【実行結果】
MCMCの実行と収束確認
MCMCを実行します。
### 事後分布からのサンプリング 10秒
with model:
idata = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
nuts_sampler='numpyro', random_seed=1234)【実行結果】省略
Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。
### r_hat>1.1の確認
# 設定
idata_in = idata # idata名
threshold = 1.1 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
収束条件を満たしています。
事後統計量を表示します。
### 推論データの要約統計情報の表示
var_names = ['mu', 'performance']
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)【実行結果】
なぜかウサギの mu の平均値はカメよりも小さくなってしまいましたｗ
モデルが良くないのでしょう。
トレースプロットを描画します。
パラメータの一部を取り上げています。
### トレースプロットの表示
pm.plot_trace(idata, var_names=var_names, compact=False)
plt.tight_layout();【実行結果】
Y の事後予測サンプリングを実行します。
### Yの事後予測サンプリングの実行
with model:
pp = pm.sample_posterior_predictive(idata, random_seed=1234)【実行結果】
ウサギの mu は峰が無く、推論を迷っている感じです。
やはり、performance の敗者＜勝者の順序を設定していないことが問題を引き起こしているようです。
これにて終了します（ジ・エンド）
10.1.5 項は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。