StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第10章「10.2.2 正の値を持つパラメータ」後編
第10章「収束しない場合の対処法」
書籍の著者 松浦健太郎 先生
この記事は、テキスト第10章「収束しない場合の対処法」・10.2節「弱情報事前分布」の 10.2.2項「正の値を持つパラメータ」の PyMC5写経 を取り扱います。
この記事は10.2.2項の後半部分、モデル式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.2.2 正の値を持つパラメータ 後編
10.2.2項の後編では、将棋のプロ棋士の勝敗データをベイズ推論します。
勝負ムラ（標準偏差）について弱情報事前分布を用いています。
インポート
### インポート
# 数値・確率計算
import pandas as pd
import numpy as np
# PyMC
import pymc as pm
import arviz as az
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')
データの読み込み・確認
サンプルコードのデータ（２ファイル）を読み込みます。
まずは勝敗データから。
### データの読み込み ◆データファイル10.3 data-shogi-player.txt
# Loser:敗者のNID、Winner:勝者のNID
data = pd.read_csv('./data/data-shogi-player.txt')
print('data.shape: ', data.shape)
display(data.head())【実行結果】
続いて棋士名データ。
### データの読み込み data-shogi-player-name.txt
data_name = pd.read_csv('./data/data-shogi-player-name.txt', delimiter='\t')
print('data_name.shape: ', data_name.shape)
display(data_name.head())【実行結果】
PyMCのモデル定義
PyMCでモデル式10-2を実装します。
モデルの定義です。
次の点でテキストとは異なるモデルになっています。
尤度関数はベルヌーイ分布。観測データは全ての値が１（勝ち）です。
「勝者パフォーマンス－敗者パフォーマンス」の値を逆ロジット関数に与えて得た「確率」を尤度関数のベルヌーイ分布に与えています。
これらをモデルに入れることで、勝者パフォーマンスが敗者パフォーマンスよりも大きくなることを企図しました。
### モデルの定義 ◆モデル式10-4 model10-4.stanモドキ
with pm.Model() as model:
### データ関連定義
## coordの定義
model.add_coord('data', values=data.index, mutable=True)
model.add_coord('player', values=data_name['kname'].values, mutable=True)
## dataの定義
# 目的変数 Y 勝者の勝ち負け区分データ（全データが勝ち:1になる）
Y = pm.ConstantData('Y', value=np.ones(len(data)), dims='data')
# 敗者インデックス 0始まりに変換
loseIdx = pm.ConstantData('loseIdx', value=data['Loser'].values - 1,
dims='data')
# 勝者インデックス 0始まりに変換
winIdx = pm.ConstantData('winIdx', value=data['Winner'].values - 1,
dims='data')
### 事前分布
# プレイヤーの強さ mu
sigmaMu = pm.Uniform('sigmaMu', lower=0, upper=100)
mu = pm.Normal('mu', mu=0, sigma=sigmaMu, dims='player')
# プレイヤーの勝負ムラ sigma
sigma = pm.Gamma('sigma', alpha=10, beta=10, dims='player')
# 各ゲームの敗者のパフォーマンス値
performanceLose = pm.Normal('performanceLose',
mu=mu[loseIdx],
sigma=sigma[loseIdx],dims='data')
# 各ゲームの勝者のパフォーマンス値
performanceWin = pm.Normal('performanceWin',
mu=mu[winIdx],
sigma=sigma[winIdx],dims='data')
# 尤度関数に用いる確率 p
p = pm.Deterministic('p', pm.invlogit(performanceWin - performanceLose),
dims='data')
### 尤度関数
obs = pm.Bernoulli('obs', p=p, observed=Y, dims='data')モデルの定義内容を見ます。
### モデルの表示
model【実行結果】
### モデルの可視化
pm.model_to_graphviz(model)【実行結果】
MCMCの実行と収束確認
MCMCを実行します。
### 事後分布からのサンプリング 4分50秒
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.09 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
収束条件を満たしています。
事後統計量を表示します。
### 推論データの要約統計情報の表示 1分10秒
var_names = ['sigmaMu', 'mu', 'sigma', 'performanceLose', 'performanceWin']
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)【実行結果】
トレースプロットを描画します。
パラメータの一部を取り上げています。
### トレースプロットの表示
var_names = ['sigmaMu', 'mu', 'sigma', 'performanceLose', 'performanceWin']
pm.plot_trace(idata, compact=False, var_names=var_names)
plt.tight_layout();【実行結果】
推定結果の解釈
強さ$${\mu}$$の推定結果について、事後平均の大きい順に一部分を表示します。
テキスト表10.1左に相当します。
### 強さμの推定結果の一部 ※事後平均meanの大きい順 ◆表10.1左
mu_stats_df = pm.summary(idata, hdi_prob=0.95, var_names=['mu'], kind='stats')
display(mu_stats_df.sort_values(by=['mean'], ascending=False).head().round(3))【実行結果】
テキストと比べて値は異なりますが、棋士はだいたい同じ結果になったと思います。
勝負ムラ$${\sigma}$$の推定結果について、事後平均の大きい順の一部分、小さい順の一部分を表示します。
テキスト表10.1右に相当します。
### 勝負ムラσの推定結果の一部 ※事後平均meanの大きい順・小さい順 ◆表10.1右
sd_stats_df = pm.summary(idata, hdi_prob=0.95, var_names=['sigma'], kind='stats')
display(sd_stats_df.sort_values(by=['mean'], ascending=False).head(3).round(3))
display(sd_stats_df.sort_values(by=['mean'], ascending=True).head(3).round(3))【実行結果】
バラツキの大きい３棋士のうち、２棋士はテキストと同じになりました。
一方で、バラツキの小さい３棋士のうち、１棋士のみテキストと同じになりました。
10.2.2 項の後編は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。