Jouez avec la modélisation statistique: quantifiez la force des équipes de la J-League avec Stan et Python

introduction

L'autre jour, j'ai lu le livre Bayes Statistical Modeling with Stan and R, qui m'a donné envie d'essayer la modélisation statistique. Je le recommande car c'est un très bon livre.

Stan peut être utilisé non seulement depuis R, mais aussi depuis la bibliothèque Python PyStan. Cette fois, je vais créer un modèle qui estime la force (comme la valeur) de chaque équipe dans la Ligue de football J1 2016 à partir des résultats du match.

La version principale de cette fois

• Python 3.5.1
• pystan==2.14.0.0

Collecter des données

Vous pouvez consulter les résultats des matchs de la J-League sur http://www.jleague.jp/stats/SFTD01.html. Supposons que vous ayez les données suivantes collectées à partir des résultats de match de chaque équipe dans la 1ère étape de 2016 d'une manière ou d'une autre. https://gist.github.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b

En pensant

• Le nombre d'équipes est de 18
• Dans la 1ère étape, chaque équipe ne jouera qu'une seule fois, donc le nombre de parties est de ** 153 ** (= 18 * 17 / 2)
• Au football, c'est presque à la maison et à l'extérieur, et il semble y avoir une différence entre le local et l'ennemi dans le match

Pensez au mécanisme par lequel les résultats de correspondance sont produits

Supposons comment le résultat du match (score) est généré comme suit.

• Supposons que chaque équipe dispose de ** puissance d'attaque ** ** puissance de défense ** ** puissance domestique **.
• Plus la puissance d'attaque est élevée, plus le score est élevé et plus la puissance de défense est élevée, plus le but est bas.
• Dans les matchs à domicile, la puissance offensive et la puissance défensive reçoivent des avantages (= puissance à domicile)

ça ira. Eh bien, cela ne semble pas cette histoire étrange pour l'instant.

Ensuite, entrons dans le monde des statistiques. Premièrement, on suppose à juste titre que le «score» suit la distribution de Poisson. En d'autres termes

Score d'une partie d'une équipe ~ Poisson (puissance d'attaque de cette équipe-puissance de défense de l'équipe adverse)

C'est comme ça. (Au fait, ~ signifie "suivre la distribution".) Cependant, les paramètres de la distribution de Poisson ne peuvent pas prendre de valeurs négatives, donc si vous prenez ʻexp` pour plus de commodité et ajoutez de l'énergie domestique,

Marquer un match pour une équipe à domicile~ Poisson(exp( 
		(Puissance domestique+Puissance d'attaque de l'équipe à domicile) -Défense de l'équipe à l'extérieur
	))

Supposer que Également symétriquement

A concédé un match dans une équipe à domicile~ Poisson(exp( 
Puissance d'attaque de l'équipe à l'extérieur- (Puissance domestique+La défense de l'équipe locale)
	))

ça ira.

Ecrire le modèle en Stan

Il s'est avéré être quelque chose comme ça.

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

explication facile

Le livre que j'ai présenté plus tôt est très utile pour Stan, mais je vais l'expliquer brièvement ici aussi.

bloc de données

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
}

Ce bloc data déclare des données externes. Les données collectées précédemment seront injectées plus tard avec ce nom.

bloc de paramètres

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

Ce bloc paramètres décrit les paramètres que vous souhaitez estimer. Cette fois, l'objectif principal est la puissance d'attaque (atk), la puissance de défense (def) et la puissance domestique (home_power). «<lower = 0>» représente la plage que le paramètre peut prendre.

bloc modèle

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]])
        ));
    }
}

Ce bloc modèle exprime la relation entre les données d'entrée et les paramètres que vous souhaitez estimer à l'aide d'une distribution de probabilité. Le "mécanisme" évoqué précédemment est également exprimé ici. De plus, on suppose que la puissance offensive et la puissance défensive de chaque équipe suivent la "distribution normale de la distribution" sigma "avec une moyenne de 0". Ce «sigma» est aussi (accessoirement) un paramètre estimé. Sinon, le calcul (MCMC) ne convergerait pas (car toutes les valeurs convergeraient dans diverses valeurs de plage). Eh bien, il semble que j'ai écrit un espoir que "je veux que vous vous mettiez dans la valeur ici pour le moment". Ce qui est intéressant, c'est qu'il était plus facile pour le calcul de converger (le «sigma» convergé était beaucoup plus petit que 1) plutôt que d'écrire le «sigma» comme un «1» fixe. Le codage en dur 1 me donne l'impression que la portée est limitée, mais il est intéressant de noter que ce n'est pas le cas.

bloc de quantités générées

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])));
        }
    }
}

Ce bloc «quantités générées» est comme un bonus. C'est une fonction qui effectue un calcul différent en utilisant la valeur convergée. Ici, les paramètres de chaque équipe sont utilisés pour échantillonner le score du match lorsque chaque équipe s'est affrontée. Vous pouvez le calculer séparément côté Python, mais comme la logique de calcul doit correspondre au modèle, il semble plus facile de le maintenir si vous le gérez ici.

Autre

Il existe un autre bloc paramètres transformés, et vous pouvez utiliser les variables data et parameters pour définir différentes variables. Des expressions telles que (home_power [th] + atk [th]) - (def [ta]) cette fois apparaissent également deux fois dans model et les quantités générées, il est donc préférable de les assembler. Je me demande, mais ...

Calculer avec le modèle et les données

Définissez quelques fonctions utiles

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'])

La fonction stan_cache met en cache le modèle stan et récupère également le code du modèle via HTTP. Pratique.

La fonction sampling_summary_to_df transforme le résultat de l'ajustement du modèle en un Pandas DataFrame. Il a l'avantage d'être plus facile à voir lorsque vous travaillez avec Jupyter.

Exécuter

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)

résultat

Un résumé des résultats de l'estimation des paramètres est par exemple le suivant. Dans le livre, "Si this Rhat est inférieur à 1,1 pour tous les paramètres estimés, on peut considérer que le calcul a convergé", donc je le considérerai comme convergé pour le moment. Cependant, il est un peu subtil que le "nombre de fois où la valeur peut être échantillonnée correctement" appelé n_eff soit inférieur à 100 (surtout lors de l'estimation de l'intervalle), donc il y a quelques subtilités, mais cette fois je le laisserai comme bon.

Voir les résultats

Estimation de la force de l'équipe

Collectez les paramètres par équipe.

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

Pour le moment, je vais les organiser par ordre décroissant de "puissance domestique".

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

est devenu. Cette fois, la force moyenne doit être de 0, donc si elle est supérieure ou inférieure à 0, c'est plus que la moyenne de l'équipe J1.

En passant, si vous le faites avec les données de la 2ème étape, vous obtiendrez un résultat complètement différent (Kashima sera bien pire). Jusqu'au dernier, c'est comme ça quand on le calcule à partir des données de la 1ère étape, donc je ne suis pas sûr.

Si vous y réfléchissez, c'est naturel, mais c'est proche du classement de la 1ère étape 2016. 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 En particulier, il semble y avoir une corrélation générale entre «le total des points et le total des buts» et «la puissance d'attaque et la puissance de défense» (bien que ce serait un problème s'il n'y en avait pas).

à la fin

J'ai essayé de simuler ce qui se passerait si j'achetais le Toto 2ème étage avec ce modèle de prédiction, mais j'ai arrêté car cela semble être un résultat terrible (rires). Il semble que nous pouvons faire de bonnes prédictions si la densité de concurrence est élevée dans un court laps de temps. Prédire 17 versets avec des données jusqu'à 16 versets, est-ce un peu mieux ??

C'est incroyable que ce genre de modélisation soit devenu facile et très amusant.