目的・モチベーション

交流戦を含むリーグ戦の試合結果をもとに、数値的に各球団の強さを表そうという試みです。

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．今後の展望

相性 (共分散を考慮する?) とか球場のアドバンテージ的なものを考慮しようかなと思っています。
その前に引き分け・・・。どうしましょう。