Stanでセ・リーグの球団の強さを推定してみた

<h2>目的・モチベーション</h2>
交流戦を含むリーグ戦の試合結果をもとに、数値的に各球団の強さを表そうという試みです。

2012年の巨人日本一を最後に、日本シリーズではパ球団が毎年優勝してますよね。正直セの方がパよりも弱いのはなんとなくお察しなのですが、どれくらい弱いんだろう。最近ベイズ統計の勉強をしているし、試しに推定してみましょうというのが、今回のモチベーションです。

注：ベイズ統計やStanの解説記事ではありません。

<h2>0．Agenda</h2>
1．モデリング
2．Stan概要
3．実装
4．推定結果
5．今後の展望

<h2>1．モデリング</h2>
試合結果から各球団の「強さ」を推定するにあたって、以下のような仮定をおきます。
まず、試合結果はより良い「パフォーマンス」を発揮した球団が勝利します。「パフォーマンス」は各球団の「強さ」に従って一定のバラつきをもって分布するランダムな要素です。

がばがばですが数式を使っていくと・・・まず、$x_t$を 球団$t$の「パフォーマンス値」とします。簡単のため、$x_t$は平均 $\mu_t$、分散 $\sigma_t^2$ の正規分布に従うものとします。$\mu_t$を球団$t$の「強さ」、$\sigma_t$を「勝負ムラ」と解釈します。
<br>
各試合の勝敗は、$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$を推定します。

<h6>(「引き分けの時どうすんの？」という疑問をお持ちの皆様。。。私もです。アドバイス頂けると幸いです。)</h6>


<h2>2．Stan概要</h2>
>Stan（スタン）はC++で書かれた統計的推論のための確率的プログラミング言語[1][2]。 Stan言語では、対数確率密度関数を計算する命令型プログラムを使用して、（ベイジアン） 統計モデルを実装できる。
<a href=https://ja.wikipedia.org/wiki/Stan_(%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%9F%E3%83%B3%E3%82%B0%E8%A8%80%E8%AA%9E)>Stan (プログラミング言語)</a>

ということでした。StanをPythonから実装する際には、<a href=https://pystan.readthedocs.io/en/latest/>Pystan</a>を使います。Pythonでベイジアンモデリングをやるなら、PyMCとかPyroとかもありますが、経済学専攻だった私にはこっちの方がしっくりきました。

<h2>3．実装</h2>
<h4>3．1 データソース</h4>
https://baseball-freak.com/game/ <br>
こちらのサイトから、スクレイピングでデータを取得しました。
こんな感じでデータを取得して・・・。

| 日付       | 球団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 | 京セラD大阪
|3月31日(金) | 日 | 西 | 1 | 8 | 18:30 | 札幌ドーム

前処理として、以下のことを行っています。
1．各球団名を1～12のインデックスに置き換え (例：巨人→1, 広島→2, など)
2．各試合ごとの勝利球団・敗北球団のインデックスを含むリストを作成

<h4>3．2 Stanコード</h4>
Stanモデルの記法はかなり直観的で、↑のモデル式をそのまま書くイメージです。

```Python
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)));
  }
}
"""
```

あとはデータを与えて事後分布のサンプリング。

```Python
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)
```

<h2> 4．結果</h2>

セ・リーグ

| 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年の交流戦順位です。

<h2> 5．今後の展望</h2>

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