統計モデリングで遊ぶ: StanとPythonでJリーグのチームの強さを数値化する

はじめに

先日 StanとRでベイズ統計モデリング という本を読みまして、統計モデリングをやってみたい気持ちになりました。とても良い本ですのでオススメです。

StanはRからだけでなく、PythonライブラリのPyStanから使うことができます。今回は、2016年のサッカーJ1リーグの各チームの強さ(のような値)を試合結果から推定するモデルを作ってみることにします。

今回の主なVersion

• Python 3.5.1
• pystan==2.14.0.0

データを集める

Jリーグの試合結果は http://www.jleague.jp/stats/SFTD01.html などから調べることができます。
そこから何らかの方法で2016年の 1stステージの 各チームの対戦結果の集めた以下のデータがあるとします。
https://gist.github.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b

away_score	away_team	date	home_score	home_team	score	visitors
1	名古屋	02/27(土)	0	磐田	0-1	14333
1	川崎Ｆ	02/27(土)	0	広島	0-1	18120
1	福岡	02/27(土)	2	鳥栖	2-1	19762
...	...	...	...	...	...	...

思うこと

• チーム数は 18 である
• 1stステージは 各チーム同士が1回だけ対戦するので試合数は 153(=18*17/2)
• サッカーではホーム＆アウェイという程だし、対戦における地元と敵地の違いがありそう

試合結果が生まれるメカニズムを考える

試合結果（スコア）がどのように生じるかを以下のように仮定してみます。

• 各チームは 攻撃力 防御力 ホーム力 というのがあるとする。
• 攻撃力が高いほど得点が多く、防御力が高いほど失点が少ない
• ホームゲームにおいては、攻撃力、防御力ともにアドバンテージが付与される(＝ホーム力)

とします。まあ、ここまではそんなにおかしな話ではなさそうです。

次に、統計の世界に入って考えてみます。
まず「得点」はポアソン分布に従うと適当に仮定します。
つまり、

あるチームのある試合の得点 ~ Poisson(そのチームの攻撃力 - 対戦チームの防御力)

という感じです。 （ちなみに ~ は「分布に従う」という意味です。）
ただし、ポアソン分布のパラメータは負の値を取れないので、expに便宜的に乗っけて、さらにホーム力を加味すると、

あるホームチームのある試合の得点 ~ Poisson(exp( 
        (ホーム力 + そのホームチームの攻撃力) - アウェイチームの防御力
    ))

になるとします。
また、対称的に

あるホームチームのある試合の失点 ~ Poisson(exp( 
        アウェイチームの攻撃力 - (ホーム力 + そのホームチームの防御力)
    ))

とします。

モデルをStanで記述する

こんな感じになりました。

model_code = """
data {
    int N;  // N Games
    int K;  // K Teams
    int Th[N]; // Home Team ID
    int Ta[N]; // Away Team ID
    int Sh[N]; // Home Team score point
    int Sa[N]; // Away Team score point
}

parameters {
    real atk[K];
    real def[K];
    real home_power[K];
    real<lower=0> sigma;
    real<lower=0> hp_sigma;
}

model {
    for (k in 1:K) {
        atk[k] ~ normal(0, sigma);
        def[k] ~ normal(0, sigma);
        home_power[k] ~ normal(0, hp_sigma);
    }

    for (n in 1:N) {
        Sh[n] ~ poisson(exp(
            (home_power[Th[n]] + atk[Th[n]]) - 
            (def[Ta[n]])
        ));

        Sa[n] ~ poisson(exp(
            (atk[Ta[n]]) - 
            (def[Th[n]] + home_power[Th[n]])
        ));
    }
}

generated quantities {
    real games[K, K, 2];

    for (th in 1:K) {
        for (ta in 1:K) {
            games[th, ta, 1] = poisson_rng(exp((home_power[th] + atk[th]) - (def[ta])));
            games[th, ta, 2] = poisson_rng(exp((atk[ta]) - (def[th] + home_power[th])));
        }
    }
}
"""

簡単な説明

Stanについては、先に紹介した本がとても参考になりますが、ここでも簡単に説明をしてみます。

dataブロック

data {
    int N;  // N Games
    int K;  // K Teams
    int Th[N]; // Home Team ID
    int Ta[N]; // Away Team ID
    int Sh[N]; // Home Team score point
    int Sa[N]; // Away Team score point
}

このdataブロックは、外部から与えられるデータを宣言しています。
先程集めたデータをこの名前で後程注入することになります。

parametersブロック

parameters {
    real atk[K];
    real def[K];
    real home_power[K];
    real<lower=0> sigma;
    real<lower=0> hp_sigma;
}

このparametersブロックは、推定したいパラメータを記述します。
今回は、攻撃力(atk)、防御力(def)、ホーム力(home_power)がメインです。
<lower=0>というのは、パラメータがとり得る範囲を表します。

modelブロック

model {
    for (k in 1:K) {
        atk[k] ~ normal(0, sigma);
        def[k] ~ normal(0, sigma);
        home_power[k] ~ normal(0, hp_sigma);
    }

    for (n in 1:N) {
        Sh[n] ~ poisson(exp(
            (home_power[Th[n]] + atk[Th[n]]) - 
            (def[Ta[n]])
        ));

        Sa[n] ~ poisson(exp(
            (atk[Ta[n]]) - 
            (def[Th[n]] + home_power[Th[n]])
        ));
    }
}

このmodelブロックは、入力データと推定したいパラメータの関係を確率分布を使って表現します。
先程の「メカニズム」もここで表現します。
また、各チームの攻撃力や防御力などは「平均0の分散sigmaの正規分布」に従うという仮定をおきます。
このsigmaも(ついでに)推定するパラメータになります。
こうしないと計算(MCMC)が収束しませんでした(全値がいろんな範囲値で収束してしまうので)。まあ「とりあえずこの辺の値で収まって欲しい」という希望を書いた感じでしょうか。
面白いのが sigma を固定の1と書くより、推定させたほうが計算が収束しやすかったことです(収束したsigmaは1よりだいぶ小さかった)。1とハードコードした方が範囲を限定した気持ちになりますが、そうでもない、というのが面白いです。

generated quantitiesブロック

generated quantities {
    real games[K, K, 2];

    for (th in 1:K) {
        for (ta in 1:K) {
            games[th, ta, 1] = poisson_rng(exp((home_power[th] + atk[th]) - (def[ta])));
            games[th, ta, 2] = poisson_rng(exp((atk[ta]) - (def[th] + home_power[th])));
        }
    }
}

このgenerated quantitiesブロックはオマケみたいなものです。
収束した値を使って、別途違う計算をしてくれる機能です。
ここでは、各チームのパラメータを使って、各チームが対戦したときの試合のスコアをサンプリングしてもらっています。別途Python側で計算してもよいのですが、計算ロジックはModelと一致させないといけないので、ここで管理したほうがメンテナンスが楽そうです。

その他

他に transformed parametersブロックというのがあって、dataやparameters の変数を使って、別変数を定義することができます。
今回の (home_power[th] + atk[th]) - (def[ta]) などの式もmodelとgenerated quantitiesに2回出てきているので、本来はまとめた方が良いのでしょうが、、、

モデルとデータで計算する

いくつか便利関数を定義しておく

from hashlib import sha256
import pickle
import time

import requests
import pandas as pd
from pystan import StanModel


def stan_cache(file=None, model_code=None, model_name=None, cache_dir="."):
    """Use just as you would `stan`"""
    # http://pystan.readthedocs.io/en/latest/avoiding_recompilation.html#automatically-reusing-models

    if file:
        if file.startswith("http"):
            model_code = requests.get(file).text
        else:
            with open(file) as f:
                model_code = f.read()
    code_hash = sha256(model_code.encode('utf8')).hexdigest()
    if model_name is None:
        cache_fn = '{}/cached-model-{}.pkl'.format(cache_dir, code_hash)
    else:
        cache_fn = '{}/cached-{}-{}.pkl'.format(cache_dir, model_name, code_hash)
    try:
        sm = pickle.load(open(cache_fn, 'rb'))
    except:
        start_time = time.time()
        sm = StanModel(model_code=model_code)
        print("compile time: %s sec" % (time.time() - start_time))
        with open(cache_fn, 'wb') as f:
            pickle.dump(sm, f)
    return sm


def sampling_summary_to_df(fit4model):
    """

    :param StanFit4model fit4model:
    :return:
    """
    s = fit4model.summary()
    return pd.DataFrame(s["summary"], columns=s['summary_colnames'], index=s['summary_rownames'])

stan_cache 関数は、stanのモデルをキャッシュしたり、モデルのコードをHTTP経由でも取得するようにしたものです。便利です。

sampling_summary_to_df関数は、modelをfitした結果を PandasのDataFrameにしてくれるものです。Jupyterなどで作業する時に見やすくなるメリットがあります。

実行する

import numpy as np
import pandas as pd

df = pd.read_csv("https://gist.githubusercontent.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b/raw/90e1c26e172b92d5f52d53d0591a611e0258bd2b/2016-j1league-1st.csv")

team_list = list(sorted(set(df['away_team']) | set(df['home_team'])))
teams = dict(zip(range(1, len(team_list)+1), team_list))
for k, v in list(teams.items()):
    teams[v] = k

model = stan_cache(model_code=model_code) 

stan_input = {
    "N": len(df),
    "K": len(team_list),
    "Th": df['home_team'].apply(lambda x: teams[x]),
    "Ta": df['away_team'].apply(lambda x: teams[x]),
    "Sh": df['home_score'],
    "Sa": df['away_score'],
}

fitted = model.sampling(data=stan_input, seed=999)
sampling_summary_to_df(fitted).sort_values("Rhat", ascending=False)

結果

パラメータの推定結果のサマリーは例えば以下のようになります。
書籍では、「このRhatが全ての推定パラメータで1.1未満になるとまあ計算が収束したとみなして良いだろう」ということですので、まあ一応収束したとみなすことにします。
ただ、n_effという「ちゃんと値がサンプリングできた回数」が100より小さいのはちょっと微妙ということ（特に区間推定するとき）なので、少し微妙なものもありますが、今回は良しとしておきます。

	mean	se_mean	sd	2.5%	25%	50%	75%	97.5%	n_eff	Rhat
lp__	-214.114891	1.542823	13.712915	-237.604900	-223.960556	-215.747425	-205.033248	-183.817981	79.0	1.049519
hp_sigma	0.094069	0.005835	0.063921	0.011408	0.044970	0.081836	0.132181	0.245524	120.0	1.031258
atk[12]	0.048430	0.009507	0.158794	-0.318161	-0.052664	0.052359	0.155878	0.342186	279.0	1.018430
atk[16]	0.225840	0.008570	0.158951	-0.075769	0.115129	0.223274	0.330571	0.531899	344.0	1.013612
sigma	0.220279	0.002456	0.056227	0.112035	0.182132	0.217550	0.255500	0.344382	524.0	1.011379
...	...	...	...	...	...	...	...	...	...	...

結果を眺める

チームの強さの推定値

チーム別にパラメータを回収します。

tk = {}

params = fitted.extract()
atks = params['atk']
defs = params['def']
home_power = params['home_power']
games = params['games']


for k, team in enumerate(team_list):
    tk[team] = {
            "atk": np.mean(atks[:, k]),
            "def": np.mean(defs[:, k]),
            "home_power": np.mean(home_power[:, k]),
    }

tdf = pd.DataFrame(tk).T

とりあえず、「ホーム力」が高い順に並べてみます。

tdf.sort_values("home_power", ascending=False)

	atk	def	home_power
鹿島	0.225840	0.217699	0.068898
浦和	0.152638	0.063192	0.041830
神戸	0.099154	-0.163383	0.030443
広島	0.293808	-0.000985	0.029861
名古屋	0.130757	-0.247969	0.015849
川崎Ｆ	0.312593	0.087859	0.014728
新潟	0.010827	-0.146252	0.011885
磐田	0.048430	-0.105230	0.011843
仙台	0.037960	-0.150151	0.005462
FC東京	-0.072115	0.017485	-0.005050
Ｇ大阪	0.087034	-0.033666	-0.005532
鳥栖	-0.232595	0.103412	-0.006186
柏	0.036172	-0.053332	-0.009568
大宮	-0.040014	0.027842	-0.020942
横浜FM	0.059420	-0.009322	-0.021393
福岡	-0.185292	-0.130759	-0.026643
甲府	0.002608	-0.268061	-0.037079
湘南	-0.000052	-0.184515	-0.049440

となりました。
今回、強さの平均が０になっているはずなので、0より高いか低いかでJ1チームの平均よりどうか、が表現されていることになります。

ちなみに、2ndステージのデータでやるとまた全然違う結果になります（鹿島がだいぶ悪くなる）。
あくまで、1stステージのデータから計算したらこうなったというだけなので、あしからず。

考えてみれば当たり前なのですが、2016年1stステージの順位表に近い感じになります。
https://ja.wikipedia.org/wiki/2016%E5%B9%B4%E3%81%AEJ1%E3%83%AA%E3%83%BC%E3%82%B0#.E9.A0.86.E4.BD.8D.E8.A1.A8
特に「総得点、総失点」と「攻撃力、防御力」は概ね相関がありそうです（無いならないで困りますが）。

さいごに

この予測モデルで2ndステージのTotoを買ったらどうなるか、、、とシミュレーションしてみようと思ったのですが、なんかひどい結果になりそうなのでやめました（笑）。もう少し、短期間で対戦密度が高ければ良い予想もできそうです。
16節までのデータで17節を予想、くらいなら少しはマシなのかな??

こういうモデリングが簡単にできるようになったというのはすごいことだし、いろいろ楽しいですね。