Skip to content

Instantly share code, notes, and snippets.

@trycycle

trycycle/npb-analysis.ipynb

Created March 16, 2020 07:34
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save trycycle/bbc9ca57143646809dcb8a7f6eb5b319 to your computer and use it in GitHub Desktop.
Save trycycle/bbc9ca57143646809dcb8a7f6eb5b319 to your computer and use it in GitHub Desktop.
Download ZIP
NumPyroによる日本プロ野球のチーム戦闘力分析
Display the source blob
Display the rendered blob
Raw
npb-analysis.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# NPBチーム戦闘力の統計モデリング\n",
"本記事では，日本プロ野球の各球団の戦闘力を統計モデリング（階層ベイズ）によって統計モデリングします．一般に，プロ野球の各球団の戦闘力は勝率，勝敗，チーム打率，チーム防御率等で表されますが，それらの指標は異なるリーグ間で比較ができません．そこで，本記事ではリーグを横断して比較できる独自指標```theta```を統計モデリングで評価します．\n",
"\n",
"今回のモデリングではPythonの確率プログラミングライブラリである[NumPyro](https://github.com/pyro-ppl/numpyro)を用います．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## データセット構築"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"本記事では，日本プロ野球機構のウェブサイトで公開されている[試合結果情報](http://npb.jp/games/)を用いて，戦闘力```theta```をモデリングします．試合結果情報ページでは，試合の結果が以下のような形式で公開されています：\n",
"\n",
"| 日時 | 球団1 | 球団1得点 | 球団2得点 | 球団2 | 球場 |\n",
"| :---: | :---: | :---: | :---: | :---: | :---: |\n",
"| 4/2 | 巨人 | 9 | 3 | 阪神 | 東京ドーム |\n",
"| 4/2 | ヤクルト | 5 | 2 | DeNA | 神宮 |\n",
"\n",
"本記事では，特定の年度（例：2019年度）の全試合結果を取得し，それらを用いてモデリングを行います．取得するデータは以下の通りです：\n",
"* 日時（year, month）\n",
"* 後攻チーム（home_team）\n",
"* 先攻チーム（road_team）\n",
"* 後攻チーム得点（home_score）\n",
"* 先攻チーム得点（road_score）\n",
"* 球場（place）\n",
"\n",
"なお，本記事ではシーズンは3月〜10月までとし，取得データはオープン戦を除く全試合（交流戦，クライマックスシリーズ，日本シリーズ含む）とします．また，分析を簡単にするために，引き分けの場合，「先攻チームの勝利」として扱います．交流戦や日本シリーズのデータを敢えて取り込んでモデリングすることで，リーグを横断して球団戦闘力を比較できるようにします．\n",
"\n",
"それでは，以下のコードを用いて，試合結果情報ページから試合結果をスクレイピングしましょう．"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# スクレイピング，データの整形に必要なライブラリ（必要に応じてインストール）\n",
"import requests\n",
"from pyquery import PyQuery as pq\n",
"import pandas as pd\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# シーズンは3月から10月までとする\n",
"REGULAR_SEASON_MONTHS = [3, 4, 5, 6, 7, 8, 9, 10] \n",
"\n",
"CENTRAL_LEAGUE = {'巨人': 'G', 'DeNA': 'DB', '阪神': 'T',\n",
" '広島': 'C', '中日': 'D', 'ヤクルト': 'S'}\n",
"\n",
"PACIFIC_LEAGUE = {'西武': 'L', 'ソフトバンク': 'H', '楽天': 'E',\n",
" 'ロッテ': 'M', '日本ハム': 'F', 'オリックス': 'Bs'}\n",
"\n",
"\n",
"def get_game_results_in_specific_month(month, year=2019):\n",
" ''' 特定年度（例: 2019年度）の特定の月（例: 5月）の全試合結果を取得する関数\n",
" '''\n",
" date = {'month':month, 'year':year}\n",
" target_url = 'http://npb.jp/games/{year}/schedule_{month:0>2}_detail.html'.format(**date)\n",
"\n",
" response = requests.get(url=target_url)\n",
" html = response.content.decode('utf8')\n",
" dom = pq(html)\n",
"\n",
" game_results = []\n",
" for tr_dom in dom(\"tr[id *= 'date']\").items():\n",
" if tr_dom('div.score1'): # 中止になっていなければ\n",
" home_team = tr_dom('div.team1').text()\n",
" road_team = tr_dom('div.team2').text()\n",
"\n",
" home_score = int(tr_dom('div.score1').text())\n",
" road_score = int(tr_dom('div.score2').text())\n",
" place = tr_dom('div.place').text().replace(\"　\", \"\")\n",
"\n",
" if home_score > road_score: # ホームチームが勝った場合\n",
" result = 1\n",
" elif home_score < road_score: # ロードチームが勝った場合\n",
" result = 0\n",
" else: # 引き分けは「ホームチームの負け」として扱う\n",
" result = 0\n",
" \n",
" # オールスター戦を除く\n",
" if home_team != 'セ・リーグ' and home_team != 'パ・リーグ':\n",
" game_results.append((year, month, home_team, road_team, home_score, road_score, result, place))\n",
" \n",
" return pd.DataFrame(game_results, columns=['year', 'month', 'home', 'road', 'home_score', 'road_score', 'result', 'place'])\n",
"\n",
"\n",
"def get_game_results_in_year(year=2019):\n",
" ''' 特定年度（例: 2019年度）の全試合結果を取得する関数．\n",
" 関数内部で上で定義したget_game_results_in_specific_monthを用いている\n",
" '''\n",
" game_result_dfs = [get_game_results_in_specific_month(month, year) for month in REGULAR_SEASON_MONTHS]\n",
" return pd.concat(game_result_dfs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"データの取得関数が用意できました．以下のコードを実行し，データフレーム```df```に2019年度の全試合の結果を格納しましょう．"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"df = get_game_results_in_year(2019)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"下記コードを実行して，取得したデータを確認してみましょう．"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>year</th>\n",
" <th>month</th>\n",
" <th>home</th>\n",
" <th>road</th>\n",
" <th>home_score</th>\n",
" <th>road_score</th>\n",
" <th>result</th>\n",
" <th>place</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>DeNA</td>\n",
" <td>中日</td>\n",
" <td>8</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>横浜</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>阪神</td>\n",
" <td>ヤクルト</td>\n",
" <td>2</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>京セラD大阪</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>広島</td>\n",
" <td>巨人</td>\n",
" <td>5</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>マツダスタジアム</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>日本ハム</td>\n",
" <td>オリックス</td>\n",
" <td>7</td>\n",
" <td>3</td>\n",
" <td>1</td>\n",
" <td>札幌ドーム</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>ロッテ</td>\n",
" <td>楽天</td>\n",
" <td>5</td>\n",
" <td>4</td>\n",
" <td>1</td>\n",
" <td>ZOZOマリン</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>ソフトバンク</td>\n",
" <td>西武</td>\n",
" <td>5</td>\n",
" <td>4</td>\n",
" <td>1</td>\n",
" <td>ヤフオクドーム</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>DeNA</td>\n",
" <td>中日</td>\n",
" <td>1</td>\n",
" <td>9</td>\n",
" <td>0</td>\n",
" <td>横浜</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>阪神</td>\n",
" <td>ヤクルト</td>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>1</td>\n",
" <td>京セラD大阪</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>広島</td>\n",
" <td>巨人</td>\n",
" <td>2</td>\n",
" <td>5</td>\n",
" <td>0</td>\n",
" <td>マツダスタジアム</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>2019</td>\n",
" <td>3</td>\n",
" <td>日本ハム</td>\n",
" <td>オリックス</td>\n",
" <td>4</td>\n",
" <td>4</td>\n",
" <td>0</td>\n",
" <td>札幌ドーム</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" year month home road home_score road_score result place\n",
"0 2019 3 DeNA 中日 8 1 1 横浜\n",
"1 2019 3 阪神 ヤクルト 2 1 1 京セラD大阪\n",
"2 2019 3 広島 巨人 5 0 1 マツダスタジアム\n",
"3 2019 3 日本ハム オリックス 7 3 1 札幌ドーム\n",
"4 2019 3 ロッテ 楽天 5 4 1 ZOZOマリン\n",
"5 2019 3 ソフトバンク 西武 5 4 1 ヤフオクドーム\n",
"6 2019 3 DeNA 中日 1 9 0 横浜\n",
"7 2019 3 阪神 ヤクルト 1 0 1 京セラD大阪\n",
"8 2019 3 広島 巨人 2 5 0 マツダスタジアム\n",
"9 2019 3 日本ハム オリックス 4 4 0 札幌ドーム"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.head(10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"結果を扱いやすくするため，得点フィールドの隣にホームチームの勝敗を示す```result```フィールドを追加しておきました．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"## 統計モデリング"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"データが得られたので，早速統計モデリングを行いましょう．\n",
"\n",
"今回のモデリングの目標は各球団の戦闘力```theta```（の分布）を推定することです．モデリングにあたって，本記事では，各試合の結果は「**対戦チームの戦闘力の差をパラメータとするベルヌーイ分布に従う」** という仮説を置きます．定式化してみましょう．今，ある試合でチーム$h$とチーム$r$が対戦したとします．チーム$h$の戦闘力，チーム$r$の戦闘力をそれぞれ $\\theta(h)$，$\\theta(r)$ とすると，ある日の試合結果 $obs$ は以下のようにモデリングできると仮定します：\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \n",
" p(h,r) = \\frac{1}{1 + e^{-(\\theta(h) - \\theta(r))}} \\\\\n",
" obs(h,r) \\sim Bernoulli(p(h,r))\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"上の式にあるように，あるチーム同士が戦った時に，ホームチーム側が勝利する確率はチームの戦闘力の差をパラメータとするシグモイド関数で定義されます．差を取っているので，チーム同士の戦闘力が拮抗しているときは，どちらが勝つかは五分五分になるようにシグモイド関数が確率を決定します．\n",
"\n",
"これでモデリングを開始，はしません．もう少し詳細な仮定を置きましょう．上の例ではチームの戦闘力は常に一定としていますが，実際には日によって戦闘力は変化します．試合に出る選手は日によって異なりますし，選手の調子も異なります．球場や天気によっても戦闘力が変わってくることも考えられます．\n",
"\n",
"そこで，本記事では，ある日のチーム $i$ の戦闘力 $\\theta(i)$ は，日によって不変のベース戦闘力 $\\theta_{base}(i)$ と試合日に依存する戦闘力 $\\theta_{delta}(i)$ から構成されると仮定します．この仮定は以下のように定式化できます：\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \n",
" \\theta_{delta}(h) \\sim Normal(0, \\sigma(h)) \\\\\n",
" \\theta(h) = \\theta_{base}(h) + \\theta_{delta}(h) \\\\\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\theta_{delta(i)}$ は平均が0，標準偏差が$\\sigma(i)$ の正規分布に従うよう仮定しています．$\\sigma(i)$ は調子のブレ具合を表現しています．これによって，試合後との調子を考慮して，その日の戦闘力をモデリングすることができます．\n",
"\n",
"\n",
"この仮定を用いて，ある年度の試合結果のモデリングを行いましょう．あるチームは$h[i]$は試合IDが$i$の試合の後攻チーム，$r[i]$は試合IDが$i$の試合の先攻チームであるとしたとき，ある年度の全試合結果を以下のようにモデリングします．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \n",
" \\theta_{base}(i) \\sim　Normal(0, 1) \\\\\n",
" \\sigma(i) \\sim HalfCauchy(1) \\\\\n",
" \\theta_{delta}(i) \\sim Normal(0, \\sigma(i)) \\\\\n",
" \\theta(i) = \\theta_{base}(i) + \\theta_{delta}(i) \\\\\n",
" p(h[i], r[i]) = \\frac{1}{1 + e^{-(\\theta(h[i]) - \\theta(r[i]))}} \\\\\n",
" obs(h[i], r[i]) \\sim Bernoulli(p(h[i], r[i]))\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 準備\n",
"では，NumPyroで上の仮定をコーディングし，MCMCサンプリングを行いましょう．まずは必要なライブラリをロードします．"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# 統計モデリング用\n",
"import numpyro\n",
"import numpyro.distributions as dist\n",
"from numpyro.diagnostics import hpdi\n",
"from numpyro.infer import MCMC, NUTS\n",
"import jax.numpy as jnp\n",
"from jax import random\n",
"\n",
"# 統計モデリングにCPU/GPUのいずれを使うかの指定\n",
"numpyro.set_platform(\"cpu\")\n",
"\n",
"# 表示桁数を小さくしておく\n",
"pd.options.display.precision = 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"分析を簡単にするために，球団名をID化しておきます．"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"team2id = {team: i for i, team in enumerate(set(df.home))}\n",
"id2team = {i: team for i, team in enumerate(set(df.home))}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"また，モデルに突っ込むためのデータを整形しておきます．"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"home_teams = jnp.array([team2id[team] for team in list(df.home)])\n",
"road_teams = jnp.array([team2id[team] for team in list(df.road)])\n",
"game_results = jnp.array(df.result)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### モデルの記述\n",
"上で述べた仮定をコーディングします．"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def model(home_teams, road_teams, results=None):\n",
" team_num = len(set(home_teams)) # チーム数\n",
" game_num = results.shape[0] # 試合数\n",
" \n",
" with numpyro.plate('team_num', team_num):\n",
" theta_base = numpyro.sample(\"theta_base\", dist.Normal(0, 1))\n",
" sigma = numpyro.sample(\"sigma\", dist.HalfCauchy(1.))\n",
"\n",
" with numpyro.plate(\"game_num\", game_num):\n",
" # 後攻チームのその日の戦闘力をサンプリング\n",
" theta_home = numpyro.sample(\"theta_home\",\n",
" dist.Normal(theta_base[home_teams], sigma[home_teams]))\n",
" \n",
" # 先攻チームのその日の戦闘力をサンプリング\n",
" theta_road = numpyro.sample(\"theta_road\",\n",
" dist.Normal(theta_base[road_teams], sigma[road_teams])) \n",
" \n",
" # 戦闘力差\n",
" p = theta_home - theta_road\n",
" \n",
" return numpyro.sample(\"obs\", dist.Bernoulli(logits=p), obs=results)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### MCMCサンプリング\n",
"設定した確率モデルに従って，MCMCによって事後分布を推定します．"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"sample: 100%|██████████| 5000/5000 [05:43<00:00, 14.55it/s, 1023 steps of size 1.68e-02. acc. prob=0.75]\n"
]
}
],
"source": [
"rng_key = random.PRNGKey(0)\n",
"rng_key, rng_key_ = random.split(rng_key)\n",
"\n",
"# バーンインは1000に指定\n",
"num_warmup, num_samples = 1000, 4000\n",
"kernel = NUTS(model)\n",
"mcmc = MCMC(kernel, num_warmup, num_samples)\n",
"\n",
"mcmc.run(rng_key, home_teams, road_teams, game_results)\n",
"\n",
"# 事後分布からのサンプル列を取得\n",
"samples = mcmc.get_samples()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 結果の表示"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"取得した事後分布のサンプルを使って，各球団のベースの戦闘力および調子のブレ（ムラッ気）を眺めてみましょう．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### ベースの戦闘力\n",
"まずはベースの戦闘力からです．事後分布からのサンプル列は変数```samples```の中に格納されています．```samples```は辞書になっており，ベースの戦闘力のサンプル列は```theta_base```をキーとして取得できます．\n",
"\n",
"以下では見やすさのために，サンプル列から計算された```theta_base```の平均値，中央値，95%HPDIをpandasのデータフレーム```theta_df```に格納しています．"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>median</th>\n",
" <th>95_low_hpdi</th>\n",
" <th>95_high_hpdi</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>ソフトバンク</th>\n",
" <td>7.058e-01</td>\n",
" <td>0.693</td>\n",
" <td>-0.072</td>\n",
" <td>1.414</td>\n",
" </tr>\n",
" <tr>\n",
" <th>西武</th>\n",
" <td>4.068e-01</td>\n",
" <td>0.403</td>\n",
" <td>-0.209</td>\n",
" <td>1.157</td>\n",
" </tr>\n",
" <tr>\n",
" <th>楽天</th>\n",
" <td>2.107e-01</td>\n",
" <td>0.210</td>\n",
" <td>-0.655</td>\n",
" <td>1.288</td>\n",
" </tr>\n",
" <tr>\n",
" <th>巨人</th>\n",
" <td>1.658e-01</td>\n",
" <td>0.143</td>\n",
" <td>-0.681</td>\n",
" <td>0.834</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ロッテ</th>\n",
" <td>1.220e-02</td>\n",
" <td>0.021</td>\n",
" <td>-1.031</td>\n",
" <td>1.023</td>\n",
" </tr>\n",
" <tr>\n",
" <th>DeNA</th>\n",
" <td>-1.275e-02</td>\n",
" <td>-0.009</td>\n",
" <td>-1.030</td>\n",
" <td>0.931</td>\n",
" </tr>\n",
" <tr>\n",
" <th>阪神</th>\n",
" <td>-2.542e-04</td>\n",
" <td>-0.030</td>\n",
" <td>-0.971</td>\n",
" <td>0.994</td>\n",
" </tr>\n",
" <tr>\n",
" <th>広島</th>\n",
" <td>-3.041e-02</td>\n",
" <td>-0.042</td>\n",
" <td>-0.923</td>\n",
" <td>0.861</td>\n",
" </tr>\n",
" <tr>\n",
" <th>日本ハム</th>\n",
" <td>-9.545e-02</td>\n",
" <td>-0.063</td>\n",
" <td>-1.041</td>\n",
" <td>0.825</td>\n",
" </tr>\n",
" <tr>\n",
" <th>オリックス</th>\n",
" <td>-2.202e-01</td>\n",
" <td>-0.239</td>\n",
" <td>-0.908</td>\n",
" <td>0.486</td>\n",
" </tr>\n",
" <tr>\n",
" <th>中日</th>\n",
" <td>-2.281e-01</td>\n",
" <td>-0.240</td>\n",
" <td>-1.020</td>\n",
" <td>0.522</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ヤクルト</th>\n",
" <td>-7.159e-01</td>\n",
" <td>-0.683</td>\n",
" <td>-1.720</td>\n",
" <td>0.205</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean median 95_low_hpdi 95_high_hpdi\n",
"ソフトバンク 7.058e-01 0.693 -0.072 1.414\n",
"西武 4.068e-01 0.403 -0.209 1.157\n",
"楽天 2.107e-01 0.210 -0.655 1.288\n",
"巨人 1.658e-01 0.143 -0.681 0.834\n",
"ロッテ 1.220e-02 0.021 -1.031 1.023\n",
"DeNA -1.275e-02 -0.009 -1.030 0.931\n",
"阪神 -2.542e-04 -0.030 -0.971 0.994\n",
"広島 -3.041e-02 -0.042 -0.923 0.861\n",
"日本ハム -9.545e-02 -0.063 -1.041 0.825\n",
"オリックス -2.202e-01 -0.239 -0.908 0.486\n",
"中日 -2.281e-01 -0.240 -1.020 0.522\n",
"ヤクルト -7.159e-01 -0.683 -1.720 0.205"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta_df = pd.DataFrame([\n",
" list(samples['theta_base'].mean(axis=0)),\n",
" list(jnp.median(samples['theta_base'], axis=0)),\n",
" list(hpdi(samples['theta_base'])[0, :]),\n",
" list(hpdi(samples['theta_base'])[1, :])\n",
"]).T.rename(columns={\n",
" 0: 'mean',\n",
" 1: 'median',\n",
" 2: '95_low_hpdi',\n",
" 3: '95_high_hpdi'\n",
"})\n",
"\n",
"theta_df.index = team2id.keys()\n",
"theta_df = theta_df.sort_values('median', ascending=False)\n",
"theta_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```theta_base```のデータフレームを眺めてみると，（ソフトバンク，西武の逆転はありますが）ベース戦闘力の中央値（median）は，セ・リーグ，パ・リーグともに大体[2019年の順位](http://npb.jp/bis/2019/stats/)に沿って並んでいます．\n",
"\n",
"今回のモデリングでは，交流戦および日本シリーズの試合結果も取り込んでモデリングを行っているため，リーグ間を横断した戦闘力が推定されています．それを踏まえて上記結果を眺めてみると，パ・リーグのAクラスチームはセ・リーグの1位（巨人）よりも戦闘力が高いとなっています．人気のセ・リーグ，実力のパ・リーグの傾向がこの分析からも見て取れます．それにしても，ヤクルトは弱いですね…\n",
"\n",
"さて，2019年のペナントレースでは，パ・リーグの勝者は西武でした．それにもかかわらず，ベースの戦闘力はソフトバンクの方が西武より上です．これはなぜでしょうか？今度は調子のブレ（$\\sigma$）を見てみましょう．以下では，サンプル列から計算された```sigma```の平均値，中央値，95%HPDIをpandasのデータフレーム```sigma_df```に格納しています．\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>mean</th>\n",
" <th>median</th>\n",
" <th>95_low_hpdi</th>\n",
" <th>95_high_hpdi</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>オリックス</th>\n",
" <td>0.615</td>\n",
" <td>0.340</td>\n",
" <td>0.039</td>\n",
" <td>1.789</td>\n",
" </tr>\n",
" <tr>\n",
" <th>中日</th>\n",
" <td>2.758</td>\n",
" <td>0.414</td>\n",
" <td>0.047</td>\n",
" <td>7.664</td>\n",
" </tr>\n",
" <tr>\n",
" <th>西武</th>\n",
" <td>0.718</td>\n",
" <td>0.518</td>\n",
" <td>0.061</td>\n",
" <td>1.431</td>\n",
" </tr>\n",
" <tr>\n",
" <th>巨人</th>\n",
" <td>1.200</td>\n",
" <td>0.575</td>\n",
" <td>0.060</td>\n",
" <td>3.414</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ソフトバンク</th>\n",
" <td>1.225</td>\n",
" <td>0.791</td>\n",
" <td>0.048</td>\n",
" <td>2.977</td>\n",
" </tr>\n",
" <tr>\n",
" <th>広島</th>\n",
" <td>6.435</td>\n",
" <td>2.478</td>\n",
" <td>0.077</td>\n",
" <td>13.620</td>\n",
" </tr>\n",
" <tr>\n",
" <th>楽天</th>\n",
" <td>6.582</td>\n",
" <td>2.746</td>\n",
" <td>0.060</td>\n",
" <td>20.699</td>\n",
" </tr>\n",
" <tr>\n",
" <th>阪神</th>\n",
" <td>5.783</td>\n",
" <td>2.839</td>\n",
" <td>0.135</td>\n",
" <td>15.267</td>\n",
" </tr>\n",
" <tr>\n",
" <th>DeNA</th>\n",
" <td>12.214</td>\n",
" <td>3.060</td>\n",
" <td>0.209</td>\n",
" <td>17.411</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ヤクルト</th>\n",
" <td>5.414</td>\n",
" <td>3.229</td>\n",
" <td>0.305</td>\n",
" <td>15.200</td>\n",
" </tr>\n",
" <tr>\n",
" <th>日本ハム</th>\n",
" <td>4.304</td>\n",
" <td>3.354</td>\n",
" <td>0.043</td>\n",
" <td>10.511</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ロッテ</th>\n",
" <td>6.425</td>\n",
" <td>4.414</td>\n",
" <td>0.149</td>\n",
" <td>14.643</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" mean median 95_low_hpdi 95_high_hpdi\n",
"オリックス 0.615 0.340 0.039 1.789\n",
"中日 2.758 0.414 0.047 7.664\n",
"西武 0.718 0.518 0.061 1.431\n",
"巨人 1.200 0.575 0.060 3.414\n",
"ソフトバンク 1.225 0.791 0.048 2.977\n",
"広島 6.435 2.478 0.077 13.620\n",
"楽天 6.582 2.746 0.060 20.699\n",
"阪神 5.783 2.839 0.135 15.267\n",
"DeNA 12.214 3.060 0.209 17.411\n",
"ヤクルト 5.414 3.229 0.305 15.200\n",
"日本ハム 4.304 3.354 0.043 10.511\n",
"ロッテ 6.425 4.414 0.149 14.643"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sigma_df = pd.DataFrame([\n",
" list(samples['sigma'].mean(axis=0)),\n",
" list(jnp.median(samples['sigma'], axis=0)),\n",
" list(hpdi(samples['sigma'])[0, :]),\n",
" list(hpdi(samples['sigma'])[1, :])\n",
"]).T.rename(columns={\n",
" 0: 'mean',\n",
" 1: 'median',\n",
" 2: '95_low_hpdi',\n",
" 3: '95_high_hpdi'\n",
"})\n",
"\n",
"sigma_df.index = team2id.keys()\n",
"sigma_df = sigma_df.sort_values('median', ascending=True)\n",
"sigma_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"中央値ベースで見ると，西武は調子のブレは3位の小ささです（オリックスはブレの小ささは1番．安定して弱い…）．一方，ソフトバンクは調子のブレでは5位．ソフトバンクは西武よりもムラッ気があると言えます．ただし，ロッテのように非常にムラッ気があるチームと比べても，調子のブレはそれほど大きくありません．むしろ戦闘力は安定している方と言えるでしょう．そこまで調子のブレも大きくなさそうですが，なぜ西武よりも戦闘力の高いソフトバンクがペナントレースで勝利数で1位になれなかったのでしょうか．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[2019年のチーム勝敗表](http://npb.jp/bis/2019/stats/std_p.html)を確認してみましたが，西武はAチームであるソフトバンク，楽天に負け越しているのに対し，ソフトバンクはAチームである西武，楽天に勝ち越しています．弱小チームから取りこぼしはあったものの（ロッテ戦ですかね？），強豪チームには勝ち越したソフトバンク．この結果が，今回の統計モデリングでは評価され，ソフトバンクの戦闘力を1位に押し上げたと考えられます．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 最後に\n",
"今回の統計モデリングでは，チームの戦闘力はベース戦闘力と試合後との調子のブレの2要素に焦点を当て，試合の勝ち負けは対戦チーム間の戦闘力の差をパラメータとするベルヌーイ分布に従うと仮定しました．その結果，上記のような結果が得られましたが，ベース戦闘力```theta_base```の推定値の分散が少し大きいのが気になりました．また，勝負には相性の要素も絡んできそうですが，今回はそれを考慮していませんでした．それらを考慮したら，もう少し厳密なモデリングができそうです（単純なモデルで現象を説明できているのであれば，それに超したことはないとは思いますが）．"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment