前回の記事で、pymc3の公式HPに、ラグビーチームの戦力をモデル化する例題があることに触れた。
https://docs.pymc.io/notebooks/rugby_analytics.html

これは、各チームの攻撃力、防御力にホーム効果の有無を付けたものになっている。これはラグビーだけでなく、サッカーなどにも使えるのでは？と思い、 Jリーグの各チームの戦力をモデルに当てはめて推定してみた。公式HPと比べた本記事の新しさは以下の2点。

J1リーグの2018年の試合を用いた。
未消化試合の結果を予測をした。

Jリーグの試合結果は、以下のサイトからダウンロードできる。 2018年のJ1全試合の結果の表をコピーし、tsvファイルとして保存した。
https://data.j-league.or.jp/SFMS01/

# データの読み込み
df = pd.read_table('./j_league/j_league_2018.tsv', delim_whitespace=False)
# columnsから空白を除去する。read_tableで、delim_whitespace=Trueにすると、未対戦結果のデータがおかしくなる。
new_col = []
for c in df.columns:
    c2 = c.strip()
    new_col.append(c2)
df.columns = new_col

df.head()

	年度	大会	節	試合日	K/O時刻	ホーム	スコア	アウェイ	スタジアム	入場者数	インターネット中継・TV放送
0	2018	Ｊ１	第１節第１日	02/23(金)	20:03	鳥栖	1-1	神戸	ベアスタ	19,633	ＤＡＺＮ
1	2018	Ｊ１	第１節第２日	02/24(土)	14:03	FC東京	1-1	浦和	味スタ	35,951	ＤＡＺＮ／ＮＨＫ　ＢＳ１／ＴＯＫＹＯ　ＭＸ（録）
2	2018	Ｊ１	第１節第２日	02/24(土)	14:03	広島	1-0	札幌	Ｅスタ	17,026	ＤＡＺＮ／ＮＨＫ広島／ＮＨＫ札幌
3	2018	Ｊ１	第１節第２日	02/24(土)	14:33	Ｇ大阪	2-3	名古屋	吹田Ｓ	28,681	ＤＡＺＮ
4	2018	Ｊ１	第１節第２日	02/24(土)	16:03	湘南	2-1	長崎	ＢＭＷス	12,148	ＤＡＺＮ

データの前処理を行う。

# ヘッダ名の変更
df = df.rename(columns={'ホーム':'home_team', 'アウェイ':'away_team'})

# チーム名をコード化する。
ch = pd.Categorical(df.home_team)
ca = pd.Categorical(df.away_team)
team_dict = dict(zip(ch.tolist(), ch.codes))
team_labes = [ x[0] for x in (sorted(team_dict.items(), key=lambda x:x[1]))]
df['i_home'] = ch.codes
df['i_away'] = ca.codes

df['スコア'].unique()

array(['1-1 ', '1-0 ', '2-3 ', '2-1 ', '0-0 ', '0-3 ', '3-3 ', '2-0 ',
       '0-1 ', '2-4 ', '2-2 ', '1-2 ', '1-3 ', '3-2 ', '3-0 ', '2-5 ',
       '3-1 ', '0-2 ', '4-4 ', '4-1 ', '4-2 ', '5-2 ', '4-0 ', '2-8 ',
       '6-2 ', '1-4 ', '中止 ', '5-1 ', 'vs '], dtype=object)

# 対戦済の結果を抽出する
df_train = df.drop(index=df[(df['スコア']=='vs ')|(df['スコア']=='中止 ')].index)

# スコアをホームチームのスコアとアウェイチームのスコアに分割する。
scores = df_train['スコア'].str.split('-', expand=True)
scores.columns = ['home_score','away_score']
df_train = pd.concat([df_train,scores], axis=1)
# 数値に変換
df_train.home_score = df_train.home_score.astype(int)
df_train.away_score = df_train.away_score.astype(int)

# 未対戦の試合を抽出する。
df_eval = df.drop(index=df[df['スコア']!='vs '].index)

対戦済結果の統計情報。やはりホームの方が有利な模様。

data = {
    'ホームが勝った確率(%)':(100*len(df_train[df_train.home_score>df_train.away_score])/len(df_train)),
    'ホームが負けた確率(%)':(100*len(df_train[df_train.home_score<df_train.away_score])/len(df_train)),
    'ホームが引分る確率(%)':(100*len(df_train[df_train.home_score==df_train.away_score])/len(df_train))
}
pd.Series(data).to_frame().T

	ホームが勝った確率(%)	ホームが負けた確率(%)	ホームが引分る確率(%)
0	42.105263	34.449761	23.444976

戦力のモデル化（公式HPと同じ）

with pm.Model() as model:
    # prior
    home = pm.Flat('home')
    intercept = pm.Flat('intercept')
    
    sd_att = pm.HalfStudentT('sd_att', nu=3, sd=2.5)
    sd_def = pm.HalfStudentT('sd_def', nu=3, sd=2.5)

    # team-specific model parameters
    atts_star = pm.Normal('atts_star', mu=0, sd=sd_att, shape=len(df.home_team.unique()))
    defs_star = pm.Normal('defs_star', mu=0, sd=sd_def, shape=len(df.home_team.unique()))

    atts = pm.Deterministic('atts', atts_star - tt.mean(atts_star))
    defs = pm.Deterministic('defs', defs_star - tt.mean(defs_star))
    
    home_theta = tt.exp(intercept+home+atts[df_train.i_home]+defs[df_train.i_away])
    away_theta = tt.exp(intercept+atts[df_train.i_away]+defs[df_train.i_home])
    
    #
    home_points = pm.Poisson('home_points', mu=home_theta, observed=df_train.home_score)
    away_points = pm.Poisson('away_points', mu=away_theta, observed=df_train.away_score)

実行時間優先のため、ADVIでパラメータ推定。

with model:
    inference = pm.ADVI()
    approx = pm.fit(n=20000, method=inference)
    trace = approx.sample(draws=5000)

Average Loss = 625.17: 100%|██████████| 20000/20000 [00:44<00:00, 452.96it/s]
Finished [100%]: Average Loss = 625.17

pm.traceplot(trace)

f:id:gaiasky:20180829221449p:plain

各チームの攻撃力と防御力。

まず、攻撃力（大きい方が点が入りやすい）。

plt.figure(figsize=(6,6))
pm.forestplot(trace, varnames=['atts'], ylabels=team_labes)

f:id:gaiasky:20180829221555p:plain

次に防御力（大きい方が点が入りやすい）。

plt.figure(figsize=(6,6))
pm.forestplot(trace, varnames=['defs'], ylabels=team_labes)

f:id:gaiasky:20180829221706p:plain

未消化試合の結果予測。

# 対戦結果のシミュレーション
def predict_vs(trace, teams,i_home, i_away):
    df_atts = pd.DataFrame(pm.stats.quantiles(trace['atts'])[50],
                           columns=['hpd_median'],
                           index=teams.i.values)
    df_defs = pd.DataFrame(pm.stats.quantiles(trace['defs'])[50],
                           columns=['hpd_median'],
                           index=teams.i.values)
    intercept = pm.stats.quantiles(trace['intercept'])[50]
    home = pm.stats.quantiles(trace['home'])[50]

    home_theta = tt.exp(intercept+home+df_atts.loc[i_home]+df_defs.loc[i_away])
    away_theta = tt.exp(intercept+df_atts.loc[i_away]+df_defs.loc[i_home])
    
    home_p = pm.Poisson.dist(mu=home_theta).random(size=1000)
    away_p = pm.Poisson.dist(mu=away_theta).random(size=1000)
    return home_p, away_p

# チーム名とindex
#以下コードは間違い。チーム名とindex番号の対応が間違っている。
#teams = df_train.home_team.unique()
#teams = pd.DataFrame(teams, columns=['team'])
#teams['i'] = teams.index
# こっちが正しい。（8/31）
teams = pd.Series(data=team_dict).to_frame(name='i').sort_values(by='i').reset_index().rename(columns={'index':'team'})

vs_pred = []
for ht,at in zip(df_eval.home_team, df_eval.away_team):
    #print('%s vs %s'%(ht,at))
    i_ht = teams[teams.team==ht].i.values[0]
    i_at = teams[teams.team==at].i.values[0]
    # シミュレーション
    hp, ap = predict_vs(trace, teams, i_ht, i_at)
    data = {
        'ホーム':ht,
        'アウェイ':at,
        'ホームの得点(期待値)':hp.mean(),
        'アウェイの得点(期待値)':ap.mean(),
        'ホームが勝つ確率(%)':100*sum(hp>ap)/len(hp),
        'ホームが負ける確率(%)':100*sum(hp<ap)/len(hp),
        'ホームが引分る確率(%)':100*sum(hp==ap)/len(hp),
    }
    vs_pred.append(pd.Series(data))
vs_pred = pd.concat(vs_pred, axis=1).T

vs_pred.head()

	ホーム	アウェイ	ホームの得点(期待値)	アウェイの得点(期待値)	ホームが勝つ確率(%)	ホームが負ける確率(%)	ホームが引分る確率(%)
0	名古屋	浦和	1.041	1.38	29.8	43.7	26.5
1	神戸	横浜FM	1.732	1.146	50.3	23.9	25.8
2	湘南	FC東京	1.135	1.243	35.2	39.6	25.2
3	鳥栖	Ｇ大阪	1.237	0.888	44.2	27	28.8
4	横浜FM	清水	1.74	1.493	43	32.9	24.1