目的・モチベーション
交流戦を含むリーグ戦の試合結果をもとに、数値的に各球団の強さを表そうという試みです。2012年の巨人日本一を最後に、日本シリーズではパ球団が毎年優勝してますよね。正直セの方がパよりも弱いのはなんとなくお察しなのですが、どれくらい弱いんだろう。最近ベイズ統計の勉強をしているし、試しに推定してみましょうというのが、今回のモチベーションです。
注：ベイズ統計やStanの解説記事ではありません。
0．Agenda
1．モデリング 2．Stan概要 3．実装 4．推定結果 5．今後の展望1．モデリング
試合結果から各球団の「強さ」を推定するにあたって、以下のような仮定をおきます。 まず、試合結果はより良い「パフォーマンス」を発揮した球団が勝利します。「パフォーマンス」は各球団の「強さ」に従って一定のバラつきをもって分布するランダムな要素です。がばがばですが数式を使っていくと・・・まず、$x_t$を 球団$t$の「パフォーマンス値」とします。簡単のため、$x_t$は平均 $\mu_t$、分散 $\sigma_t^2$ の正規分布に従うものとします。$\mu_t$を球団$t$の「強さ」、$\sigma_t$を「勝負ムラ」と解釈します。
各試合の勝敗は、$x_{t1} > x_{t2}$ の時$t1$ 球団の勝利となり、それ以外で$t2$ の勝利です。$x_{t1}$ と $x_{t2}$が独立に分布する (相性とか想定しない) 場合、$x_{t1}-x_{t2} ～ Normal(\mu_{t1}-\mu_{t2},\sigma_{t1}^2 + \sigma_{t2}^2)$ なので、$t1$ が$t2$ に勝利する確率は以下の式で表されます $P(x_{t1} - x_{t2} > 0 | \mu_{t1}, \mu_{t2}, \sigma_{t1}, \sigma_{t2}) $。
以下では実際のデータを元にパラメータ$\mu$や$\sigma$を推定します。
(「引き分けの時どうすんの？」という疑問をお持ちの皆様。。。私もです。アドバイス頂けると幸いです。)
2．Stan概要
>Stan（スタン）はC++で書かれた統計的推論のための確率的プログラミング言語[1][2]。 Stan言語では、対数確率密度関数を計算する命令型プログラムを使用して、（ベイジアン） 統計モデルを実装できる。 Stan (プログラミング言語)ということでした。StanをPythonから実装する際には、Pystanを使います。Pythonでベイジアンモデリングをやるなら、PyMCとかPyroとかもありますが、経済学専攻だった私にはこっちの方がしっくりきました。
3．実装
3．1 データソース
https://baseball-freak.com/game/こちらのサイトから、スクレイピングでデータを取得しました。 こんな感じでデータを取得して・・・。
| 日付 | 球団1 | 球団2 | 球団1得点 | 球団2得点 | 球場 |
|---|---|---|---|---|---|
| 3月31日(金) | 広 | 阪 | 6 | 10 | 18:00 |
| 3月31日(金) | 巨 | 中 | 6 | 2 | 18:00 |
| 3月31日(金) | ヤ | D | 9 | 2 | 18:00 |
| 3月31日(金) | オ | 楽 | 4 | 6 | 18:30 |
| 3月31日(金) | 日 | 西 | 1 | 8 | 18:30 |
前処理として、以下のことを行っています。
1．各球団名を1～12のインデックスに置き換え (例：巨人→1, 広島→2, など)
2．各試合ごとの勝利球団・敗北球団のインデックスを含むリストを作成
3．2 Stanコード
Stanモデルの記法はかなり直観的で、↑のモデル式をそのまま書くイメージです。stan_code = """
data {
int G; // 総試合数
int T; // 総球団数 = 12
int winner_idx[G]; // 勝利球団のインデックス
int loser_idx[G]; // 敗北球団のインデックス
}
parameters {
vector[T] mu_pef; // 各球団のパフォーマンス平均 = 強さ
vector<lower=0>[T] s_pef; // 各球団のパフォーマンスのばらつき = 勝負ムラ
}
model {
mu_pef ~ normal(0, 1); // パフォーマンス平均の事前分布
s_pef ~ gamma(10, 10); // パフォーマンスばらつきの事前分布
for (i in 1:G) {
target += normal_lcdf(0 |
mu_pef[loser_idx[i]] - mu_pef[winner_idx[i]],
sqrt(pow(s_pef[winner_idx[i]], 2) + pow(s_pef[loser_idx[i]], 2)));
}
}
"""
あとはデータを与えて事後分布のサンプリング。
import pystan
data = {
'G': len(loser_idx),
'T': len(teams_idx),
'winner_idx': winner_idx,
'loser_idx': loser_idx,
}
sm = pystan.StanModel(model_code=stan_code)
fit = sm.sampling(data=data, iter=1000, chains=2, warmup=200, thin=1, seed=101,n_jobs=-1)
4．結果
セ・リーグ
| 2019年順位 | 球団名 | 強さ | 勝負ムラ |
|---|---|---|---|
| 1位 | 巨人 | 0.062 | 1.428 |
| 2位 | DeNA | -0.235 | 1.345 |
| 3位 | 阪神 | -0.252 | 1.301 |
| 4位 | 広島 | -0.266 | 1.362 |
| 5位 | 中日 | -0.317 | 1.1403 |
| 6位 | ヤクルト | -0.588 | 1.374 |
パ・リーグ
| 2019年順位 | 球団名 | 強さ | 勝負ムラ |
|---|---|---|---|
| 1位 | 西武 | 0.165 | 1.235 |
| 2位 | ソフトバンク | 0.122 | 1.143 |
| 3位 | 楽天 | -0.028 | 1.455 |
| 4位 | ロッテ | -0.076 | 1.495 |
| 5位 | 日ハム | -0.172 | 1.223 |
| 6位 | オリックス | -0.271 | 1.125 |
(ちゃんと収束してんの？というツッコミはとりあえずおいておいてください)
基本的にセ・パ共にリーグ戦順位の順にパフォーマンス平均 (強さ) が高いということがわかりました。
なんとなく直観的な結果な気がします。
セ・リーグでは、リーグ優勝した巨人の値が頭一つ抜けていますね。2位のDeNAともゲーム差がそれなりにあったので、なんとなく納得です。
それにしてもヤクルト・・・ずば抜けて・・・。セ・パの比較をしてみると、セ・リーグ優勝の巨人も、パ・リーグの中では2位のソフトバンクより弱いという結果でした。日本シリーズでも負けちゃったしなぁ。
| 2019年交流戦順位 | 球団名 | 強さ |
|:-----------|------------:|:------------:|:------------:|
| 1位 | ソフトバンク | 0.122 | 1.052 |
| 2位 | オリックス | -0.271 | 1.055 |
| 3位 | 巨人 | 0.062 | 1.122 |
| 4位 | DeNA | -0.235 | 1.159 |
| 5位 | 西武 | 0.165 | 1.062 |
| 6位 | 楽天 | -0.028 | 1.026 |
| 7位 | 日ハム | -0.172 | 1.052 |
| 8位 | 中日 | -0.317| 1.055 |
| 9位 | ロッテ | -0.076 | 1.122 |
| 10位 | 阪神 | -0.266 | 1.159 |
| 11位 | ヤクルト | -0.588 | 1.062 |
| 12位 | 広島 | -0.266 | 1.026 |
参考までに、↑が2019年の交流戦順位です。
5．今後の展望
相性 (共分散を考慮する?) とか球場のアドバンテージ的なものを考慮しようかなと思っています。
その前に引き分け・・・。どうしましょう。