第６章「男子プロテニスプレイヤーの強さモデリング」のベイズモデリングをPyMC Ver.5 で
この記事は、テキスト「たのしいベイズモデリング２」の第６章「男子プロテニスプレイヤーの強さモデリング」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。
この章では、男子プロテニスプレイヤーの勝敗データを用いて「強さ」を推論します。
テキストでは２つのモデルを論じています。
そして今回も、自己流PyMCにとって実装の困難度が高いモデル。
芳しい結果を出すことはできませんでした（泣）
まあまあ堅苦しいことはさておき、PyMCのベイズモデリングの世界を楽しみましょう。
テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。
テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。
サマリー
テキストの概要
執筆者 ： 伊藤瑛志 先生
モデル難易度： ★★★★★ （難しい）
自己評価
評点
$$
\begin{array}{c:c:c}
実装精度 & ★★・・・ & いまいち \\
結果再現度 & ★・・・・ & NG!! \\
楽しさ & ★★★★・ & 楽しそう \\
\end{array}
$$
評価ポイント
最大の難関はStanの「Orderd型」の実装です。
この型は「敗者パフォーマンス＜勝者パフォーマンス」の大小関係を確保できるものです。
PyMCでの実装方法が分からず（注）、仕方ないので、モデルを改造しました。
１つのモデルは収束しましたが、もう１つは収束しませんでした。
（注）
PyMCのtransforms.orderedによる順序制約は１次元配列には効くものの、２次元配列には効かないような雰囲気です。
工夫・喜び・反省
なんと！テキストのサンプルコードはPythonなのです！型破り！（注）
おかげさまで、円滑にコードを読み取ることができました。
なお、サンプルコードの一部に理解できない実装箇所があり、結果として、テキストと異なるデータ処理を採用しました。
（注）まえがきをご覧いただき、型破りをご堪能ください。
北大路書房のWebサイト（URL省略）の本書サポートページから入手できるデータとフリーの統計解析環境RとStanのコードによって、本書の内容はすべて再現できる。
モデルの概要
テキストの調査・実験の概要
■ プロテニスプレイヤーの勝敗結果のデータセット
ネットで公開されているATP（男子プロテニス協会）ツアーの勝敗結果データセットを利用します。
データセットを提供する Jeff Sackmann 氏に感謝です。
この記事でもJeff氏のATPデータセットを引用いたします。
データソースは次のとおりです。
■ 分析対象プレイヤーと対象期間
テキストは２つのモデルを構築しています。
最初のモデルでは2015年1月～2017年2月の勝敗データを用いて、次の18プレイヤーを分析します。
R・フェデラー
R・ナダル
N・ジョコビッチ
A・マレー
S・ワウリンカ
J・M・デル=ポトロ
M・ラオニッチ
錦織圭
G・モンフィス
T・ベルディヒ
JW・ツォンガ
D・フェレール
R・ガスケ
M・チリッチ
G・ディミトロフ
D・ティーム
N・キリオス
A・ズベレフ
最初のモデルでは2005年1月～2016年12月の勝敗データを用いて、次の9プレイヤーを分析します。
R・フェデラー
R・ナダル
N・ジョコビッチ
A・マレー
S・ワウリンカ
J・M・デル=ポトロ
錦織圭
T・ベルディヒ
D・フェレール
数年前のデータですが、超有名プレイヤーが名を連ねていて、結果が楽しみですね！
【お知らせ】
研究・調査を実施していないようですので、今回は図示を省略いたします。
テキストのモデリング
テキストは２つのモデルを提案しています。
１．強さのモデル
やっぱり図がみたいですよね！（作りたいだけ）
この「１．強さのモデル」のイメージを図示しました。
図を描いたら気持ちがスッキリしました。
モデルの話を続けます。
■目的変数と関心のあるパラメータ
目的変数は・・・なんと！不明です！
目的変数の観測値を制約にもつ尤度関数はありません。
関心のあるパラメータは、プレイヤー$${n}$$の強さ$${\mu_n}$$と勝負ムラ$${\sigma_{pf\ n}}$$です。
$$
\begin{align*}
performance_{g,1} &\sim \text{Normal}\ (\mu_{loser},\ \sigma_{pf\ loser}), \quad g=1, \cdots, G \\
performance_{g,2} &\sim \text{Normal}\ (\mu_{winner}\, \sigma_{pf\ winner}), \quad g=1, \cdots, G \\
performance_{g,1} &< performance_{g,2}, \quad g=1, \cdots, G \\
\mu_n &\sim \text{Normal}\ (0,\ \sigma{\mu}), \quad n=1, \cdots, N \\
\sigma_{pf\ n} &\sim \text{Gamma}\ (10,\ 10), \quad n=1, \cdots, N \\
\end{align*}
$$
$${performance}$$は１回のゲーム$${g}$$で発揮する力です。
$${performance_{g,1}}$$は敗者$${loser}$$のパフォーマンス、$${performance_{g,2}}$$は勝者$${winner}$$のパフォーマンスです。
$${performance_{g,1} < performance_{g,2}}$$は、あるゲームの敗者のパフォーマンスが勝者のパフォーマンスよりも小さいという制約条件です。
Stanではorderd型を用いることで、このパフォーマンス間の制約条件を実現できます。PyMCの実装はと言うと・・・波乱含みです！
２．時系列で見た強さのモデル
時系列で各プレイヤーの強さがどのように推移したのかを見るモデルです。
目的変数、関心のあるパラメータは１つ目のモデルと同じです。
$$
\begin{align*}
performance_{g,1y} &\sim \text{Normal}\ (\mu_{loser_y},\ \sigma_{pf\ loser_y}), \quad g=1, \cdots, G, \quad y=1, \cdots, Y \\
performance_{g,2y} &\sim \text{Normal}\ (\mu_{winner_y}\, \sigma_{pf\ winner_y}), \quad g=1, \cdots, G, \quad y=1, \cdots, Y \\
performance_{g,1y} &< performance_{g,2y}, \quad g=1, \cdots, G, \quad y=1, \cdots, Y \\
\mu_{n_1} &\sim \text{Normal}\ (0,\ \sigma_{\mu_{n_1}}), \quad n=1, \cdots, N \\
\mu_{n_y} &\sim \text{Normal}\ (\mu_{n_{y-1}},\ \sigma_{\mu_{n_{y-1}}}), \quad n=1, \cdots, N, \quad y=2, \cdots, Y \\
\sigma_{pf\ n_y} &\sim \text{Gamma}\ (10,\ 10), \quad n=1, \cdots, N, \quad y=1, \cdots, Y \\
\sigma_{\mu_{n_y}} &\sim \text{Normal}\ (0,\ 1), \quad n=1, \cdots, N, \quad y=1, \cdots, Y \\
\end{align*}
$$
時間軸$${y}$$が追加されました。
また、$${\sigma_{\mu_{n_y}}}$$の従う分布が明記されました。
テキストのモデルは以上となります。
ひとつお知らせがあります。
自己流PyMCモデルでは、テキストの忠実な再現ができませんでした。
やむを得ず、改造モデルを構築することにしました。
PyMC実装の章に改造の詳細を書きます。
■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
PyMCの自己流モデルはテキストとパラメータ値の尺度が変わったり、収束しなかったりと、様々な要因があり、推論結果で分析するには厳しい状況です。
PyMC実装
Let's enjoy PyMC & Python !
準備・データ確認
１．インポート
### インポート
# ユーティリティ
import pickle
# 数値・確率計算
import pandas as pd
import numpy as np
# PyMC
import pymc as pm
import arviz as az
# WEBアクセス
import requests
# データ前処理
from sklearn.preprocessing import LabelEncoder
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')２．データの読み込みと前処理
JeffさんのWebサイトからCSVファイルをダウンロードして分析に用います。
まずダウンロードしましょう。
「path1」に「CSVファイルの格納フォルダ」を指定します。
ファイルは年単位に作成されています。
### データのダウンロード 25秒
## 初期値設定
# 取得する年の範囲（2005年から2017年まで）
years = list(range(2005, 2018))
# ダウンロードファイルの格納フォルダ
path1 = r'./data/'
# ダウンロードファイルのURL（ファイル名を除く）
url1 = r'https://github.com/JeffSackmann/tennis_atp/blob/master/'
# 2005～2017年のcsvファイルを１つずつダウンロードして格納フォルダに保存
for year in years:
# csvファイル名の設定
filename = rf'atp_matches_{year}.csv'
# ファイル取得URLの設定
url = url1 + filename + r'?raw=true'
# ダウンロード格納パスの設定
path = path1 + filename
# Web接続
res = requests.get(url)
# ダウンロードファイルの保存
with open(path, mode='wb') as f:
f.write(res.content)【実行結果】なし
続いて、CSVファイルをpandasのデータフレームに読み込みます。
サンプルコード preprocess.py を参考にして、一部コードを書き換えて作成しました。
年単位のファイルをマージして、分析期間の絞り込みを行っています。
### データの作成 preprocess.pyを参考にして書き換え
## csvファイルを読み込んで結合
df = pd.concat([pd.read_csv(path1 + f'atp_matches_{year}.csv',
usecols=[5, 10, 18]) for year in years])
## データの前処理１
# tourney_dateに欠損値があるレコードを削除（実際には欠損値は無い）
df = df.dropna(subset=['tourney_date'])
# 「年」列の追加
df['year'] = df['tourney_date'].apply(lambda x: int(str(x)[:4]))
# date列をdatetime型に変換
df['date'] = pd.to_datetime(df['tourney_date'].apply(lambda x: str(x)[:8]))
# データの期間を2017年2月までに絞り込み ★テキストは未代入のため絞り込みできていない
df = df[df['date'] < '2017-03-01']
# dfの列を限定してインデックスを再設定
df = df[['year', 'winner_name', 'loser_name']].reset_index(drop=True)
## データの表示
display(df)
## データの保存
df.to_csv('data.csv', index=False)【実行結果】
ゲームの年 year、勝者名 winner_name、敗者名 loser_name です。
次の箇所はテキストのサンプルコードを採用しないで、改造しています。
df = df[df['date'] < '2017-03-01']テキストの場合、「df =」の記載がないので、もしかすると、2017年3月1日以降のゲームを含んで分析しているかもしれません。
モデル構築１「強さのモデル」
モデルごとにデータ作成とモデル構築を実践します。
１．データの前処理
サンプルコード preprocess.py を参考にして、一部コードを書き換えて作成しました。
### 前処理２ analyze_strength.pyを参考に書き換え
## 初期値設定
# プレイヤーリスト（英語） ★Wawrinka, del Potro, Tsongaの表記が変更されている
arr_target_player = np.array([
'Roger Federer', 'Rafael Nadal', 'Novak Djokovic', 'Andy Murray',
'Stan Wawrinka', 'Juan Martin del Potro', 'Milos Raonic',
'Kei Nishikori', 'Gael Monfils', 'Tomas Berdych', 'Jo-Wilfried Tsonga',
'David Ferrer', 'Richard Gasquet', 'Marin Cilic', 'Grigor Dimitrov',
'Dominic Thiem', 'Nick Kyrgios', 'Alexander Zverev'
])
# プレイヤーリスト（日本語）
arr_target_player_ja = np.array([
'R・フェデラー', 'R・ナダル', 'N・ジョコビッチ', 'A・マレー',
'S・ワウリンカ', 'J・M・デル=ポトロ', 'M・ラオニッチ',
'錦織圭', 'G・モンフィス', 'T・ベルディヒ', 'JW・ツォンガ',
'D・フェレール', 'R・ガスケ', 'M・チリッチ', 'G・ディミトロフ',
'D・ティーム', 'N・キリオス', 'A・ズベレフ'
])
## dfから2015年以降、注目１８選手を抽出してdataを作成
# 抽出条件の設定 ★テキストは二重に抽出している
query = ((df['year'] >= 2015)
& (df['winner_name'].isin(arr_target_player))
& (df['loser_name'].isin(arr_target_player)))
# データの抽出
data = df[query]
## LabelEncoderを利用してプレイヤーIDを算出し、勝者・敗者にセット
# ラベルエンコーダーのクラスにarr_target_playerを設定
# arr_target_playerの順番、0始まりでIDを生成
le = LabelEncoder()
le.classes_ = arr_target_player
# 勝者名（英語）にラベルエンコーダーを適用して勝者のプレイヤーIDを設定
data['winner_id'] = le.transform(data['winner_name'])
# 敗者名（英語）にラベルエンコーダーを適用して敗者のプレイヤーIDを設定
data['loser_id'] = le.transform(data['loser_name'])
## 日本語のプレイヤー名をセット
data['勝者'] = data['winner_id'].apply(lambda x: arr_target_player_ja[x])
data['敗者'] = data['loser_id'].apply(lambda x: arr_target_player_ja[x])
## 列順の変更、インデックスの再取得
data = data[['year', 'winner_name', 'loser_name', '勝者', '敗者',
'winner_id', 'loser_id']].reset_index(drop=True)
## 前処理後データの表示
display(data)【実行結果】
368ゲームのデータになりました。
次の２箇所はテキストのサンプルコードを採用しないで、改造しています。
まず１箇所目。
# プレイヤーリスト（英語） ★Wawrinka, del Potro, Tsongaの表記が変更されている
arr_target_player = np.array([
'Roger Federer', 'Rafael Nadal', 'Novak Djokovic', 'Andy Murray',
'Stan Wawrinka', 'Juan Martin del Potro', 'Milos Raonic',
'Kei Nishikori', 'Gael Monfils', 'Tomas Berdych', 'Jo-Wilfried Tsonga',
'David Ferrer', 'Richard Gasquet', 'Marin Cilic', 'Grigor Dimitrov',
'Dominic Thiem', 'Nick Kyrgios', 'Alexander Zverev'
])プレイヤーの英名がテキスト執筆時点から変わっている感じです。
Wawrinka選手、del Potro選手、Tsonga選手の英名を変更しました。
もう１箇所はこちら。
## dfから2015年以降、注目１８選手を抽出してdataを作成
# 抽出条件の設定 ★テキストは二重に抽出している
query = ((df['year'] >= 2015)
& (df['winner_name'].isin(arr_target_player))
& (df['loser_name'].isin(arr_target_player)))
# データの抽出
data = df[query]
テキストの抽出条件は個人的に使いづらそうだったので、こちらの条件に変更しました。
２．データの外観の確認
各プレイヤーの勝ち数、負け数を棒グラフで見てみましょう。
棒グラフの描画には、seabornのbarplotを利用しています。
### 勝数、負数の棒グラフの描画
## 棒グラフ用のデータ加工
# 勝者の勝ち数の算出
win = data['勝者'].value_counts()
# 敗者の負け数の算出
lose = data['敗者'].value_counts()
## 描画
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
# 勝者の棒グラフの描画
sns.barplot(x=win.keys(), y=win.values, color='tab:blue', ax=ax1)
ax1.tick_params(axis='x', labelrotation=70)
ax1.set_ylabel('勝ち・負け回数')
# 敗者の棒グラフの描画
sns.barplot(x=lose.keys(), y=lose.values, color='tomato', ax=ax2)
ax2.tick_params(axis='x', labelrotation=70)
plt.tight_layout();【実行結果】
準決勝・決勝などの上位のゲームに進んだプレイヤーは、戦ったゲーム数が多くて、しかも、ゲームに勝っているからゲーム数が多いという状況です。
３．モデルの再考
テキストのモデルをそのまま実装することができなかったので（個人のスキルの問題）、モデル自体を大きく改造しました。
実装の障害になったのは「StanのOrderd型」です。
ゲームごとに「敗者パフォーマンス＜勝者パフォーマンス」の大小関係を定義することは、私の力量では困難でした。
ここで、自己流モデルの説明をいたします。
まず、尤度関数は「勝者が勝つ＝値１」を目的変数・観測値とする「ベルヌーイ分布」です。
ベルヌーイ分布の確率パラメータには、「勝者パフォーマンス－敗者パフォーマンス」の値をロジスティック関数（逆ロジット関数、シグモイド関数）で確率値に変換した値を用いました。
この２本を用いることで、「勝者パフォーマンス－敗者パフォーマンス」が正の値（勝者パフォーマンスの方が大きい値）になることを狙っています。
数式にすると、こんな感じ。
観測値$${obs}$$はすべての値が$${1}$$です。
$$
\begin{align*}
obs &\sim \text{Bernoulli}\ (p) \\
p &= \text{logistic}\ (performance_{winner}-performance_{loser}) \\
\end{align*}
$$
もう１つのモデルもベルヌーイ分布＋ロジステッィク関数を用いています。
モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。
$$
\begin{align*}
\sigma_{\mu} &\sim \text{HalfNormal}\ (\text{sigma}=1) \\
\mu &\sim \text{HalfNormal}\ (\text{sigma}=\sigma_{\mu}) \\
\sigma &\sim \text{Gamma}\ (\text{alpha}=12,\ \text{beta}=8) \\
performance_{loser} &\sim \text{Normal}\ (\text{mu}=\mu_{loser},\ \text{sigma}=\sigma) \\
performance_{winner} &\sim \text{Normal}\ (\text{mu}=\mu_{winner},\ \text{sigma}=\sigma) \\
p &= \text{invlogit}\ (performance_{winner}-performance_{loser}) \\
likelihood &\sim \text{Bernoulli}\ (\text{p}=p) \\
\end{align*}
$$
もう１点、テキストと変えている箇所があります。
勝負ムラ$${\sigma}$$が従うガンマ分布のパラメータ$${\alpha,\beta}$$の値です。
$${\sigma}$$は正規分布の標準偏差パラメータに用いるので「正の値」にする必要があります。
ガンマ分布は確率変数が正の値をとるので、標準偏差にも使えるようです。
標準偏差の分布によく使う半正規分布、半コーシー分布と比べて、パラメータ値のありかたを検討します。
まずはガンマ分布の確率密度関数を可視化します。
### 様々なパラメータのガンマ分布の描画
# パラメータ値のリスト化
alphabeta = [(3, 3), (4, 4), (5, 5), (8,8), (10, 10), (12, 12), (20, 20), (12, 8)]
# x軸の値
testx = np.linspace(0, 4, 1001)
# 様々なパラメータのガンマ分布を描画
for (alpha, beta) in alphabeta:
# ガンマ分布の確率密度関数を算出 scipy.statsのgamma
testy = stats.gamma.pdf(testx, a=alpha, scale=1/beta)
# ガンマ分布の確率密度関数の描画
plt.plot(testx, testy, label=rf'$\alpha$={alpha}, $\beta$={beta}')
plt.legend();【実行結果】
パラメータ$${\alpha, \beta}$$に同一値を与えて比較します。
青色がパラメータ値３。右に裾の長い分布です。
パラメータ値を大きくするにつれて、峰の値が大きくなり、かつ左右対称に近づいています。
今回採用した$${\alpha=12, \beta=8}$$（グレーの曲線）は、尖りを抑えて幅広な値を許容するものになります。
テキストのパラメータ値は$${\alpha=10, \beta=10}$$（紫の曲線）です。
半正規分布、半コーシー分布を見てみましょう。
### 半正規分布、半コーシー分布の例
## 描画データの作成
# x軸の値
testx = np.linspace(0, 4, 1001)
# 半正規分布の確率密度関数
test_halfnorm = stats.halfnorm.pdf(testx, loc=0.5, scale=1)
# 半コーシー分布の確率密度関数
test_halfcauchy = stats.halfcauchy.pdf(testx, loc=0.5, scale=1)
## 描画
# 半正規分布の確率密度関数の描画
plt.plot(testx, test_halfnorm, label=rf'HalfNormal(loc=0.5, $\sigma$=1)')
# 半コーシー分布の確率密度関数の描画
plt.plot(testx, test_halfcauchy, label=rf'HalfCauchy(loc=0.5, $\gamma$=1)')
plt.legend();【実行結果】
曲線左側が切り立っています。
ガンマ分布のほうが緩やかな曲線を描いていて、峰の左右の分布を幅広くすくい上げられるような感じです。
４．モデルの定義
数式表現を実直にモデル記述します。
### モデルの定義
with pm.Model() as model:
### coordの定義
model.add_coord('game', values=data.index, mutable=True)
model.add_coord('player', values=arr_target_player_ja, mutable=True)
model.add_coord('losewin', values=['lose', 'win'], mutable=True)
### dataの定義
# 観測値を仮作成：勝者の勝敗フラグ（すべて1：勝ち）
y = pm.ConstantData('y', value=np.ones(len(data)), dims='game')
# 敗者インデックス
loseIdx = pm.ConstantData('loseIdx', value=data['loser_id'].values,
dims='game')
# 勝者インデックス
winIdx = pm.ConstantData('winIdx', value=data['winner_id'].values,
dims='game')
### 事前分布
# 強さμ
sigmaMu = pm.HalfNormal('sigmaMu', sigma=1)
mu = pm.HalfNormal('mu', sigma=sigmaMu, dims='player')
# 勝負ムラσ
# sigma = pm.Gamma('sigma', alpha=10, beta=10, dims='player')
sigma = pm.Gamma('sigma', alpha=12, beta=8, dims='player')
### パフォーマンス値
# 敗者のパフォーマンス値
performanceLose = pm.Normal('performanceLose',
mu=mu[loseIdx],
sigma=sigma[loseIdx],
dims='game')
# 勝者のパフォーマンス値
performanceWin = pm.Normal('performanceWin',
mu=mu[winIdx],
sigma=sigma[winIdx],
dims='game')
### 勝者の成功確率
# （勝者パフォーマンス － 敗者パフォーマンス）のロジスティック関数
p = pm.Deterministic('p', pm.invlogit(performanceWin - performanceLose),
dims='game')
### 尤度 観測値はすべて1（全ゲームが成功1）
likelihood = pm.Bernoulli('likelihood', p=p, observed=y, dims='game')【モデル注釈】
coordの定義
座標に名前を付けたり、その座標が取りうる値を設定できます。
今回は次の３つを設定しました。データの行の座標：名前「game」、値「行インデックス」
プレイヤーの座標：名前「player」、値「プレイヤーの日本語名」
配列の座標：名前「losewin」、値「lose, win」
dataの定義
次の３つを設定しました。目的変数・観測値に用いるすべての値が１の観測値 y
敗者インデックス loseIdx
勝者インデックス winIdx
パラメータの事前分布
$${\sigma_{\mu}}$$は区間広めの一様分布を採用しました。
その他のパラメータはテキストの記載のとおりです。
尤度
観測値 y は確率 p をパラメータにするベルヌーイ分布に従います。
５．モデルの外観の確認
# モデルの表示
model【実行結果】
いろんな分布が使われています。
# モデルの可視化
pm.model_to_graphviz(model)【実行結果】
左の角丸四角の中は対決感があります！
６．事後分布からのサンプリング
乱数生成数はテキストと同じです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ40秒でした。
### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用 40秒
# テキスト：iter=1000, warmup=指定なし(たぶん500), chains=4,
with model:
idata = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
nuts_sampler='numpyro', random_seed=2019)【実行結果】
７．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$とします。
### r_hat>1.1の確認
# 設定
idata_in = idata # idata名
threshold = 1.05 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
$${\hat{R} > 1.05}$$でカウントしたところ、０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。
ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。
### 推論データの要約統計情報の表示
pm.summary(idata, hdi_prob=0.95).round(3)【実行結果】
続いてトレースプロットを確認しましょう。
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['mu'])
plt.tight_layout();【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['sigma'])
plt.tight_layout();【実行結果】
収束している感じです。
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['sigmaMu'])
plt.tight_layout();【実行結果】
収束している感じです。
分析 強さモデル
関心のあるパラメータ「強さ $${\mu_n}$$」と「勝負ムラ $${\sigma_{pf\ n}}$$」をテキストにならって分析します。
①強さ $${\mu_n}$$
テキストの図6.1で示される「強さ$${\mu_n}$$」をプロットします。
seaborn の violinplot を利用します。
### μの事後分布のヴァイオリンプロットの描画 ※図6.1に相当
# 事後分布サンプリングデータからμを取り出し
mu_data = idata.posterior.mu.stack(sample=('chain', 'draw')).data.T
# 描画
plt.figure(figsize=(10, 4))
# ヴァイオリンプロットの描画
sns.violinplot(mu_data)
# 修飾：x軸目盛ラベルの回転
plt.xticks(range(0, 18), arr_target_player_ja, rotation=60);【実行結果】
テキストよりも全体的に値が大きくなりました。
（テキストは０～４の範囲、こちらは０～８）
ジョコビッチ選手がダントツで強く、また、左の４プレイヤーは他のプレイヤーよりも強さが引き立っています。
また、テキストよりもデル＝ポトロ選手の強さが強調されている感じになりました。
続いて、テキスト表6.2の強さ$${\mu_n}$$の事後中央値、不偏標準偏差、50%CI区間を算出します。
まず事後統計量の計算関数を定義します。
### 事後分布統計量の計算関数の定義
# 中央値、不偏標準偏差、50%信用区間を算出
def calc_stats(x):
return np.median(x), np.std(x, ddof=1), np.quantile(x, 0.25), \
np.quantile(x, 0.75)【実行結果】なし
続いて事後統計量を計算します。
### μの事後分布統計量の計算、データフレーム化
# データフレームの初期化
mu_stats_df = pd.DataFrame()
# 18人のプレイヤーごとに統計量計算を繰り返し処理
for player in arr_target_player_ja:
# 指定プレイヤーのμの値を平坦化して取得
tmp_data = idata.posterior.mu.sel(player=player).data.flatten()
# 指定プレイヤーのμの事後分布統計量を計算
tmp_stats = calc_stats(tmp_data)
# データフレームに追加
mu_stats_df = pd.concat([mu_stats_df, pd.DataFrame(tmp_stats).T], axis=0)
# データフレームの最終化（インデックス名、カラム名の設定）
mu_stats_df.index = arr_target_player_ja
mu_stats_df.columns = ['Median', 'STD', '25%', '75%']
# データフレームの表示
display(mu_stats_df.round(3))【実行結果】
テキストよりも不偏標準偏差の値が大きいです。
ばらつきが大きくなりました。
ではテキスト表6.2のようなヒートマップ形式の表示をしてみます。
seaborn の heatmap を利用します。
### μの事後分布統計量の可視化 ※表6.2に相当
# ４つの統計量ごとにヒートマップの描画を繰り返し処理
for i, col in enumerate(mu_stats_df.columns):
# １列目のヒートマップにだけプレイヤー名を表示
yticklabels = False
if i == 0:
yticklabels = True
# axesの確保
ax = plt.subplot(1, 4, i+1)
# ヒートマップの描画
sns.heatmap(data=mu_stats_df[[col]], ax=ax, annot=True, fmt='.3f',
cmap='Oranges', cbar=False, yticklabels=yticklabels)
# 統計量名を上に表示
ax.xaxis.tick_top()
plt.tight_layout();【実行結果】
テキストの中央値では、強い順にジョコビッチ選手、フェデラー選手、マレー選手ですが、こちらはジョコビッチ選手、マレー選手、フェデラー選手の順になりました。
続いて、テキストで「確認のため」に用いられる表6.3の勝敗戦績ヒートマップを描画します。
まずデータを加工してpandasのデータフレームにまとめます。
行のプレイヤーが列のプレイヤーに何勝したかを集計します。
### 対戦組み合わせごとの勝敗数をデータフレーム化
# NumPy配列を初期化（18×18、0埋め）
game_mtx = np.zeros((len(arr_target_player_ja), len(arr_target_player_ja)),
dtype=np.int32)
# dataデータフレームの勝者IDと敗者IDの組み合わせごとに勝敗数をカウント
for (i, j) in zip(data['winner_id'], data['loser_id']):
game_mtx[i, j] += 1
# データフレーム化
game_mtx_df = pd.DataFrame(game_mtx, index=arr_target_player_ja,
columns=arr_target_player_ja)
# データフレームの表示
display(game_mtx_df)【実行結果】
描画します！
### 2015年1月～2017年2月の各プレイヤーの勝敗戦績のヒートマップを描画 ※表6.3に相当
# ★テキストは2017年4月～12月のゲームを含んでいるため値が不一致となる
# ヒートマップの描画（上のセルで作成したデータフレームを利用）
plt.figure(figsize=(7, 7))
ax = sns.heatmap(data=game_mtx_df, annot=True, cmap='Reds', cbar=False,
linewidths=0.5)
# 修飾：x軸目盛ラベルの回転
ax.tick_params(axis='x', labelrotation=70)
# 修飾：x軸ラベル、y軸ラベルの表示
ax.set(xlabel='負けプレイヤー', ylabel='勝ちプレイヤー');【実行結果】
この図の「勝ち数の大きさ」（色が濃い）を眺めてみると、フェデラー選手よりもマレー選手の方がよく勝っている感じがします。
なお、このヒートマップの数値はテキストと異なっています。
理由は、データの前処理時に分析対象期間がテキストと異なっているからです。
②勝負ムラ $${\sigma_{pf\ n}}$$
テキストの図6.2で示される「勝負ムラ$${\sigma_{pf\ n}}$$」をプロットします。
seaborn の violinplot を利用します。
### σの事後分布のヴァイオリンプロットの描画 ※図6.2に相当
# 事後分布サンプリングデータからσを取り出し
sigma_data = idata.posterior.sigma.stack(sample=('chain', 'draw')).data.T
# 描画
plt.figure(figsize=(10, 4))
# ヴァイオリンプロットの描画
sns.violinplot(sigma_data)
# 修飾：x軸目盛ラベルの回転
plt.xticks(range(0, 18), arr_target_player_ja, rotation=60);【実行結果】
どのプレイヤーも似たようなばらつきになり、分析に使いにくいものになってしまいました。
強さ値の広がりが$${\mu_n}$$に反映されてしまったのかもしれません。
続いて、強さのときと同様に勝負ムラの事後統計量をヒートマップで表示します。
まず描画に利用するデータを作成してデータフレーム化します。
### σの事後分布統計量の計算、データフレーム化
# データフレームの初期化
sigma_stats_df = pd.DataFrame()
# 18人のプレイヤーごとに統計量計算を繰り返し処理
for player in arr_target_player_ja:
# 指定プレイヤーのσの値を平坦化して取得
tmp_data = idata.posterior.sigma.sel(player=player).data.flatten()
# 指定プレイヤーのσの事後分布統計量を計算
tmp_stats = calc_stats(tmp_data)
# データフレームに追加
sigma_stats_df = pd.concat([sigma_stats_df, pd.DataFrame(tmp_stats).T],
axis=0)
# データフレームの最終化（インデックス名、カラム名の設定）
sigma_stats_df.index = arr_target_player_ja
sigma_stats_df.columns = ['Median', 'STD', '25%', '75%']
# データフレームの表示
display(sigma_stats_df.round(3))【実行結果】
中央値 Median はどんぐりの背比べになっています。
描画します。
### σの事後分布統計量の可視化 ※表6.4に相当
# ４つの統計量ごとにヒートマップの描画を繰り返し処理
for i, col in enumerate(mu_stats_df.columns):
# １列目のヒートマップにだけプレイヤー名を表示
yticklabels = False
if i == 0:
yticklabels = True
# axesの確保
ax = plt.subplot(1, 4, i+1)
# ヒートマップの描画
sns.heatmap(data=sigma_stats_df[[col]], ax=ax, annot=True, fmt='.3f',
cmap='Oranges', cbar=False, yticklabels=yticklabels)
# 統計量名を上に表示
ax.xaxis.tick_top()
plt.tight_layout();【実行結果】
値の大小関係はテキストに近い感じがします。
マレー選手は各ゲームで勝ち負けのムラが一番少ないようです。
最後に、勝利確率$${p}$$の事後サンプリングデータの状況を確認しましょう。
目論見的には $${p>0.5}$$を狙っていました。
### 勝利確率pの事後分布サンプリングデータのヒストグラムの描画
plt.hist(idata.posterior.p.stack(sample=('chain', 'draw')).mean(axis=1).data,
bins=15)
plt.axvline(0.5, color='red');【実行結果】
大半は 0.5 を超えています。
しかし、0.5 未満のケースも若干あります。
print('0.5以下の件数: ',
((idata.posterior.p.stack(sample=('chain', 'draw'))
.mean(axis=1).data) <= 0.5).sum())【実行結果】
0.5 未満のケースは 8 件でした。
この 8 件が強さや勝負ムラに大きく影響しているとは考えにくいです。
テキストの結果とこちらの結果の違いは、直感的には、ロジスティック関数を用いて「勝者パフォーマンス」と「敗者パフォーマンス」の差を決めてしまう過程が要因のように感じます。
モデル構築２「時系列で見た強さのモデル」
１．データの前処理
サンプルコード preprocess.py を参考にして、一部コードを書き換えて作成しました。
### データ前処理 analyze_strength_ts.pyを参考に書き換え
### csvファイルの読み込み
data2 = pd.read_csv('data.csv')
## 初期値設定
# プレイヤーリスト（英語）★Wawrinka, del Potroの表記が変更されている
arr_target_player = np.array([
'Roger Federer', 'Rafael Nadal', 'Novak Djokovic',
'Andy Murray', 'Stan Wawrinka', 'Juan Martin del Potro',
'Kei Nishikori', 'Tomas Berdych', 'David Ferrer'
])
# プレイヤーリスト（日本語）
arr_target_player_ja = np.array([
'R・フェデラー', 'R・ナダル', 'N・ジョコビッチ',
'A・マレー', 'S・ワウリンカ', 'J・M・デル=ポトロ',
'錦織圭', 'T・ベルディヒ', 'D・フェレール'
])
# 年のリスト
years2 = list(range(2005, 2017))
## 2005年～2016年の注目９選手を抽出してdataを作成
# 抽出条件の設定 ★（データ検証していないが）テキストは二重に抽出していると思われる
query = ((data2['year'] <= 2016)
& (data2['winner_name'].isin(arr_target_player))
& (data2['loser_name'].isin(arr_target_player)))
# データの抽出
data2 = data2[query]
## LabelEncoderを利用して年ID（0始まり）を算出
# ラベルエンコーダーのクラスにyearを設定
le = LabelEncoder()
le.classes_ = np.array(years2)
# year列にラベルエンコーダーを適用して年IDを設定
data2['year_id'] = le.transform(data2['year'])
## LabelEncoderを利用してプレイヤーIDを算出し、勝者・敗者にセット
# ラベルエンコーダーのクラスにarr_target_playerを設定
# arr_target_playerの順番、0始まりでIDを生成
le = LabelEncoder()
le.classes_ = arr_target_player
# 勝者名（英語）にラベルエンコーダーを適用して勝者のプレイヤーIDを設定
data2['winner_id'] = le.transform(data2['winner_name'])
# 敗者名（英語）にラベルエンコーダーを適用して敗者のプレイヤーIDを設定
data2['loser_id'] = le.transform(data2['loser_name'])
## 日本語のプレイヤー名をセット
data2['勝者'] = data2['winner_id'].apply(lambda x: arr_target_player_ja[x])
data2['敗者'] = data2['loser_id'].apply(lambda x: arr_target_player_ja[x])
## 列順の変更、インデックスの再取得
data2 = data2[['year', 'year_id', 'winner_name', 'loser_name', '勝者', '敗者',
'winner_id', 'loser_id']].reset_index(drop=True)
## 前処理後データの表示
display(data2)【実行結果】
660ゲームのデータになりました。
２．データの外観の確認
各プレイヤーの勝ち数、負け数を棒グラフで見てみましょう。
棒グラフの描画には、seabornのbarplotを利用しています。
### 勝数、負数の棒グラフの描画
## 棒グラフ用のデータ加工
# 勝者の勝ち数の算出
win = data2['勝者'].value_counts()
# 敗者の負け数の算出
lose = data2['敗者'].value_counts()
## 描画
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
# 勝者の棒グラフの描画
sns.barplot(x=win.keys(), y=win.values, color='tab:blue', ax=ax1)
ax1.tick_params(axis='x', labelrotation=70)
ax1.set_ylabel('勝ち・負け回数')
# 敗者の棒グラフの描画
sns.barplot(x=lose.keys(), y=lose.values, color='tomato', ax=ax2)
ax2.tick_params(axis='x', labelrotation=70)
plt.tight_layout();【実行結果】
３．モデルの再考
「強さモデル」と同様にベルヌーイ分布＋ロジスティック関数を用いたモデルに改造しました。
モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。
$$
\begin{align*}
\sigma_{\mu} &\sim \text{HalfNormal}\ (\text{sigma}=1) \\
init\_dist &\sim \text{HalfNormal}\ (\text{sigma}=\sigma_{\mu}) \\
\mu &\sim \text{GaussianRandomWalk}\ (\text{sigma}=\sigma_{\mu},\ \text{init\_dist}=init\_dist) \\
\sigma &\sim \text{Gamma}\ (\text{alpha}=12,\ \text{beta}=8) \\
performance_{loser} &\sim \text{Normal}\ (\text{mu}=\mu_{loser},\ \text{sigma}=\sigma) \\
performance_{winner} &\sim \text{Normal}\ (\text{mu}=\mu_{winner},\ \text{sigma}=\sigma) \\
p &= \text{invlogit}\ (performance_{winner}-performance_{loser}) \\
likelihood &\sim \text{Bernoulli}\ (\text{p}=p) \\
\end{align*}
$$
４．モデルの定義
数式表現を実直にモデル記述します。
強さ$${\mu_{n_y}}$$の時系列には GaussianRandomWalk を用いました。
### モデルの定義
with pm.Model() as model2:
### coordの定義
model2.add_coord('game', values=data2.index, mutable=True)
model2.add_coord('year', values=years2, mutable=True)
model2.add_coord('player', values=arr_target_player_ja, mutable=True)
### dataの定義
# 観測値を仮作成：勝者の勝敗フラグ（すべて1：勝ち）
y = pm.ConstantData('y', value=np.ones(len(data2)), dims='game')
# 年インデックス
yearIdx = pm.ConstantData('yearIdx', value=data2['year_id'].values,
dims='game')
# 敗者インデックス
loseIdx = pm.ConstantData('loseIdx', value=data2['loser_id'].values,
dims='game')
# 勝者インデックス
winIdx = pm.ConstantData('winIdx', value=data2['winner_id'].values,
dims='game')
### 事前分布
# μの正規分布の標準偏差 ※dimsが他のパラメータと逆になる
sigmaMu = pm.HalfNormal('sigmaMu', sigma=1, dims=('year', 'player'))
# μの0期の初期分布
init_dist = pm.HalfNormal.dist(sigma=sigmaMu[0, :])
# 強さμ：performanceが従う正規分布の平均
mu = pm.GaussianRandomWalk('mu', mu=0, sigma=sigmaMu[1:, :],
init_dist=init_dist, dims=('player', 'year'))
# 強さの変化具合σ：performanceが従う正規分布の標準偏差
sigma = pm.Gamma('sigma', alpha=12, beta=8, dims=('player', 'year'))
### パフォーマンス値
# 敗者のパフォーマンス値
performanceLose = pm.Normal('performanceLose',
mu=mu[loseIdx, yearIdx],
sigma=sigma[loseIdx, yearIdx],
dims='game')
# 勝者のパフォーマンス値
performanceWin = pm.Normal('performanceWin',
mu=mu[winIdx, yearIdx],
sigma=sigma[winIdx, yearIdx],
dims='game')
### 勝者の成功確率
# （勝者パフォーマンス － 敗者パフォーマンス）のロジスティック関数
p = pm.Deterministic('p', pm.invlogit(performanceWin - performanceLose),
dims='game')
### 尤度 観測値はすべて1（全ゲームについて成功1を設定）
likelihood = pm.Bernoulli('likelihood', p=p, observed=y, dims='game')【モデル注釈】
coordの定義
座標に名前を付けたり、その座標が取りうる値を設定できます。
今回は次の３つを設定しました。データの行の座標：名前「game」、値「行インデックス」
ゲーム年の座標：名前「yaer」、値「ゲーム開催年」
プレイヤーの座標：名前「player」、値「プレイヤーの日本語名」
dataの定義
次の３つを設定しました。目的変数・観測値に用いるすべての値が１の観測値 y
敗者インデックス loseIdx
勝者インデックス winIdx
パラメータの事前分布
$${\mu_{n_y}}$$はGaussianRandomWalkに従います。
その他のパラメータはテキストの記載のとおりです。
尤度
観測値 y は確率 p をパラメータにするベルヌーイ分布に従います。
５．モデルの外観の確認
# モデルの表示
model2【実行結果】
いろんな分布が使われています。
# モデルの可視化
pm.model_to_graphviz(model2)【実行結果】
変数の次元にゲーム開催年 year=12 が現れています。
６．事後分布からのサンプリング
乱数生成数はテキストと同じです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ5分50秒でした。
### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用 5分50秒
# テキスト：iter=5000, warmup=指定なし(たぶん2500), chains=4
with model2:
idata2 = pm.sample(draws=2500, tune=2500, chains=4, target_accept=0.99,
nuts_sampler='numpyro', random_seed=2019)【実行結果】
７．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$とします。
### r_hat>1.1の確認
# 設定
idata_in = idata2 # idata名
threshold = 1.1 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
$${\hat{R} > 1.1}$$のパラメータが多数存在します。
収束条件を満たすことができませんでした。
収束できていない推論データでは分析に使えません（泣）
そこで、ここからは、一般的な推論データ確認と、分析には使えませんがテキストのコードをPython化して動かしてみます。
手始めにざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。
### 推論データの要約統計情報の表示
pm.summary(idata2, hdi_prob=0.95, var_names=['mu', 'sigma', 'sigmaMu'])【実行結果】
続いてトレースプロット（一部のパラメータ）を確認しましょう。
### トレースプロットの表示
pm.plot_trace(idata2, compact=False, var_names=['mu'])
plt.tight_layout();【実行結果】
左のグラフの４本のマルコフ連鎖chainはそれぞれ異なる軌跡を描いています。
右のグラフでは線が満遍なく描画されておりません。
グラフ下部は真っ黒なバーコード（発散の記号）で埋め尽くされています。
収束していません。
分析に至らないデータ確認 時系列で見た強さモデル
テキストにならって、MCMCサンプリングデータを可視化していきましょう！
収束していないので、データの中身は適切ではありません。
ご留意くださいませ。
まずはテキストの図6.3「ジョコビッチ選手と錦織圭選手の時系列プロット」です。
### ジョコビッチ選手と錦織選手の強さμの事後分布の中央値・50%HDI区間の描画
# ※図6.3に相当
# ジョコビッチ選手と錦織選手の年別中央値の算出
Djoko = (idata2.posterior.mu.stack(sample=('chain', 'draw'))
.sel(player=['N・ジョコビッチ']).squeeze().median(axis=1))
Nishi = (idata2.posterior.mu.stack(sample=('chain', 'draw'))
.sel(player=['錦織圭']).squeeze().median(axis=1))
# ジョコビッチ選手と錦織選手の50%HDI区間の算出
Djoko_hdi = (az.hdi(idata2.posterior.mu
.sel(player=['N・ジョコビッチ']), hdi_prob=0.5)
.mu.squeeze().data)
Nishi_hdi = (az.hdi(idata2.posterior.mu
.sel(player=['錦織圭']), hdi_prob=0.5)
.mu.squeeze().data)
# 描画
plt.figure(figsize=(10, 5))
# ジョコビッチ選手の中央値の描画
plt.plot(years2, Djoko, label='N・ジョコビッチ選手')
# ジョコビッチ選手の50$HDI区間の描画
plt.fill_between(years2, Djoko_hdi[:, 0], Djoko_hdi[:, 1], alpha=0.2)
# 錦織選手の中央値の描画
plt.plot(years2, Nishi, label='錦織圭選手')
# 錦織選手の50%HDI区間の描画
plt.fill_between(years2, Nishi_hdi[:, 0], Nishi_hdi[:, 1], alpha=0.2)
# 修飾
plt.xlabel('年')
plt.ylabel('強さ')
plt.legend(loc='upper left')
plt.show()【実行結果】
線は事後中央値、塗りつぶしは50%CI区間です。
錦織圭選手の中央値が下方にブレてしまいました。
テキストはマイナス値になっていません。
続いてテキストの図6.4「強さ$${\mu_{n_y}}$$の事後中央値の時系列プロット」です。
### μの事後分布の中央値の時系列プロットの描画 ※図6.4に相当
# 事後分布サンプリングデータから9選手の中央値を抽出
median9 = idata2.posterior.mu.stack(sample=('chain', 'draw')).median(axis=2)
# 描画
plt.figure(figsize=(10, 5))
# プレイヤーごとに描画処理
for player in arr_target_player_ja:
# 各プレイヤーの時系列折れ線グラフの描画
plt.plot(years2, median9.sel(player=player).data, label=f'{player}')
# 凡例をグラフの外に表示
plt.legend(bbox_to_anchor=(1, 1))
# 修飾
plt.xlabel('年')
plt.ylabel('強さ')
plt.show()【実行結果】
上位選手はテキストと「なんとなく」近い感じがします。
続いて、テキストの図6.5、勝負ムラ$${\sigma_{pf\ n_y}}$$の描画です。
描画結果は全くテキストと異なります。
### σの事後分布のヴァイオリンプロットの描画 ※図6.5に相当
# 事後分布サンプリングデータからσを取り出し
mu_data2 = idata2.posterior.sigma.stack(sample=('chain', 'draw'))
# 描画
fig, ax = plt.subplots(3, 3, figsize=(10, 10), sharey=True)
for i, player in enumerate(arr_target_player_ja):
# axesの設定
axes = divmod(i, 3)
# ヴァイオリンプロットの描画
tmp = mu_data2.sel(player=player).data.T
sns.violinplot(data=tmp, ax=ax[axes])
# 修飾：x軸目盛ラベルの回転
ax[axes].set_xticks(range(12), years2, rotation=60)
# 修飾：タイトル（プレイヤー名）の表示
ax[axes].set_title(player)
plt.tight_layout()
plt.show()【実行結果】
最後のプロットは、テキストの図6.6、「強さの変化具合$${\sigma_{\mu_{n_y}}}$$の事後分布のプロット」です。
この描画結果も全くテキストと異なります。
### σ_μの事後分布のヴァイオリンプロットの描画 ※図6.6に相当
# 事後分布サンプリングデータからσ_μを取り出し
sigma_data2 = idata2.posterior.sigmaMu.stack(sample=('chain', 'draw'))
# 描画
fig, ax = plt.subplots(3, 3, figsize=(10, 10), sharey=True)
for i, player in enumerate(arr_target_player_ja):
# axesの設定
axes = divmod(i, 3)
# ヴァイオリンプロットの描画
tmp = sigma_data2.sel(player=player).data.T
sns.violinplot(data=tmp, ax=ax[axes])
# 修飾：x軸目盛ラベルの回転
ax[axes].set_xticks(range(12), years2, rotation=60)
# 修飾：タイトル（プレイヤー名）の表示
ax[axes].set_title(player)
plt.tight_layout()
plt.show()【実行結果】
以上で分析もどき「データのプロット」を終了します。
最後に、推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
idataをpickleで保存します。
### idataの保存 pickle
file = r'idata_ch06.pkl'
with open(file, 'wb') as f:
pickle.dump(idata, f)
file = r'idata2_ch06.pkl'
with open(file, 'wb') as f:
pickle.dump(idata2, f)読み込みコードは次のとおりです。
### idataの読み込み pickle
file = r'idata_ch06.pkl'
with open(file, 'rb') as f:
idata_load = pickle.load(f)
### idataの読み込み pickle
file = r'idata2_ch06.pkl'
with open(file, 'rb') as f:
idata2_load = pickle.load(f)以上で第６章は終了です。
おわりに
記憶に残る章
テキスト「たのしいベイズモデリング２」のページをペラペラめくっていたとき、最初に惹かれたのがこの第６章でした。
私も名前を知っているメジャーなテニスプレイヤーが取り扱われていたこと、そしてゲーム勝敗からプレイの強さを推論できることに惹かれました。
そこで、テキストで最初に第６章のPyMC化を実践したのです。
いざモデリングをしてみると、めっちゃクセの強いモデルであること、めっちゃPyMC化しにくいことが分かりました。
目的変数が観測値じゃないぞ！
StanのOrderd型ってどうすればいいの！？
Pythonのスクリプトファイル（.pyファイル）を用いているとは言え、処理途中経過の出力が無いので、コード内で何やっているのか良く分からん・・・
特に、StanのOrderd型の代替を探すために、PyMC公式サイトのさまざまなコードを眺めては自環境で動かしてみる、といった試行錯誤を繰り返し実践したので、時間がとてもかかってしまいました。
しかし、最適な代替策を見つけることができませんでした（泣）
難しかったゆえに記憶に残る章になりました。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？