Spielen Sie mit statistischer Modellierung: Quantifizieren Sie die Stärke von J-League-Teams mit Stan und Python

Einführung

Neulich las ich das Buch Bayes Statistical Modeling mit Stan und R, das mich dazu brachte, statistische Modellierung auszuprobieren. Ich empfehle es, weil es ein sehr gutes Buch ist.

Stan kann nicht nur von R aus verwendet werden, sondern auch von der Python-Bibliothek PyStan. Dieses Mal werde ich ein Modell erstellen, das die Stärke (den gleichen Wert) jeder Mannschaft in der Soccer J1 League 2016 anhand der Spielergebnisse schätzt.

Diesmal ist die Hauptversion

• Python 3.5.1
• pystan==2.14.0.0

Daten sammeln

Sie können die Ergebnisse der J-League-Spiele unter http://www.jleague.jp/stats/SFTD01.html überprüfen. Angenommen, Sie haben die folgenden Daten von jedem Team in der ersten Phase des Jahres 2016 in irgendeiner Weise gesammelt. https://gist.github.com/mokemokechicken/a26eb7230e51fba7018476f7706e5b0b

Denken

• Die Anzahl der Teams beträgt 18
• In der 1. Phase wird jedes Team nur einmal spielen, daher beträgt die Anzahl der Spiele ** 153 ** (= 18 * 17/2)
• Im Fußball geht es um Heim und Auswärts, und es scheint einen Unterschied zwischen dem lokalen und dem feindlichen Land im Spiel zu geben.

Denken Sie über den Mechanismus nach, mit dem Übereinstimmungsergebnisse erzeugt werden

Nehmen wir an, wie das Match-Ergebnis (Score) wie folgt generiert wird.

• Angenommen, jedes Team hat ** Offensivkraft ** ** Defensivkraft ** ** Heimkraft **.
• Je höher die Angriffskraft, desto höher die Punktzahl und je höher die Verteidigungskraft, desto niedriger das Ziel.
• In Heimspielen erhalten sowohl Offensivkraft als auch Defensivkraft Vorteile (= Heimkraft)

Wird besorgt. Nun, es scheint bisher nicht so seltsam zu sein.

Als nächstes betreten wir die Welt der Statistik. Zunächst wird angemessen angenommen, dass die "Punktzahl" der Poisson-Verteilung folgt. Mit anderen Worten

Punktzahl eines Spiels eines Teams ~ Poisson (Angriffskraft dieses Teams - Verteidigungskraft des gegnerischen Teams)

Es ist wie es ist. (Übrigens bedeutet "~" "der Verteilung folgen".) Die Parameter der Poisson-Verteilung können jedoch keine negativen Werte annehmen. Wenn Sie also der Einfachheit halber die Ausgangsleistung zu "exp" hinzufügen,

Ein Spiel für eine Heimmannschaft erzielen~ Poisson(exp( 
		(Heimstrom+Angriffskraft der Heimmannschaft) -Auswärtsmannschaft Verteidigung
	))

Nehme an, dass Auch symmetrisch

Ein Spiel in einer Heimmannschaft kassiert~ Poisson(exp( 
Angriffskraft des Auswärtsteams- (Heimstrom+Die Verteidigung der Heimmannschaft)
	))

Wird besorgt.

Schreiben Sie das Modell in Stan

Es stellte sich heraus, dass es so etwas war.

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

einfache Erklärung

Für Stan ist das Buch, das ich zuvor vorgestellt habe, sehr hilfreich, aber ich werde es auch hier kurz erklären.

Datenblock

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
}

Dieser "Daten" -Block deklariert extern gegebene Daten. Die früher gesammelten Daten werden später mit diesem Namen versehen.

Parameterblock

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

Dieser Parameterblock beschreibt die Parameter, die Sie schätzen möchten. Dieses Mal liegt der Schwerpunkt auf Angriffskraft (atk), Verteidigungskraft (def) und Heimkraft (home_power). <lower = 0> repräsentiert den Bereich, den der Parameter einnehmen kann.

Modellblock

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

Dieser "Modell" -Block drückt die Beziehung zwischen den Eingabedaten und den Parametern aus, die Sie mithilfe einer Wahrscheinlichkeitsverteilung schätzen möchten. Der zuvor erwähnte "Mechanismus" wird auch hier ausgedrückt. Außerdem wird angenommen, dass die Offensiv- und Defensivkraft jedes Teams der "Normalverteilung des Verteilungssigmas" mit einem Durchschnitt von 0 "folgt. Dieses "Sigma" ist auch (im Übrigen) ein geschätzter Parameter. Andernfalls würde die Berechnung (MCMC) nicht konvergieren (da alle Werte in verschiedenen Bereichswerten konvergieren würden). Nun, es scheint, als hätte ich eine Hoffnung geschrieben, dass "ich möchte, dass Sie vorerst in den Wert hier passen". Interessant ist, dass es für die Berechnung einfacher war, zu konvergieren (das konvergierte "Sigma" war viel kleiner als 1), als das "Sigma" als feste "1" zu schreiben. Durch die Hardcodierung 1 habe ich das Gefühl, dass der Bereich begrenzt ist, aber es ist interessant, dass dies nicht der Fall ist.

generierte Mengen blockieren

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

Dieser Block "generierte Mengen" ist wie ein Bonus. Diese Funktion berechnet anhand des konvergierten Werts anders. Hier werden die Parameter jeder Mannschaft verwendet, um die Punktzahl des Spiels abzutasten, wenn jede Mannschaft gegeneinander spielt. Sie können es auf der Python-Seite separat berechnen. Da die Berechnungslogik jedoch mit dem Modell übereinstimmen muss, scheint es einfacher zu sein, es zu verwalten, wenn Sie es hier verwalten.

Andere

Es gibt einen weiteren Block "transformierte Parameter", und Sie können die Variablen "Daten" und "Parameter" verwenden, um verschiedene Variablen zu definieren. Ausdrücke wie "(home_power [th] + atk [th]) - (def [ta])" erscheinen dieses Mal auch zweimal in "model" und "generierten Mengen", daher ist es besser, sie zusammenzusetzen. Ich frage mich, aber ...

Berechnen Sie mit Modell und Daten

Definieren Sie einige nützliche Funktionen

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

Die Funktion "stan_cache" speichert das Stan-Modell zwischen und ruft den Modellcode auch über HTTP ab. Praktisch.

Die Funktion sample_summary_to_df wandelt das Ergebnis der Anpassung des Modells in einen Pandas DataFrame um. Es hat den Vorteil, dass es bei der Arbeit mit Jupyter leichter zu erkennen ist.

Lauf

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)

Ergebnis

Eine Zusammenfassung der Parameterschätzungsergebnisse lautet beispielsweise wie folgt. In dem Buch "Wenn dieses" Rhat "für alle geschätzten Parameter kleiner als 1,1 ist, kann davon ausgegangen werden, dass die Berechnung konvergiert hat", daher werde ich davon ausgehen, dass es vorerst konvergiert. Es ist jedoch ein wenig subtil, dass die "Häufigkeit, mit der der Wert richtig abgetastet werden kann" mit der Bezeichnung "n_eff" weniger als 100 beträgt (insbesondere beim Schätzen des Intervalls), daher gibt es einige Feinheiten, aber dieses Mal werde ich es als gut belassen.

Zeigen Sie die Ergebnisse an

Geschätzte Teamstärke

Sammeln Sie die Parameter nach Team.

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

Vorerst werde ich sie in absteigender Reihenfolge der "Heimatmacht" anordnen.

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

ist geworden. Dieses Mal sollte die durchschnittliche Stärke 0 sein. Ob sie also höher oder niedriger als 0 ist, ist mehr als der Durchschnitt des J1-Teams.

Übrigens, wenn Sie es mit den Daten der 2. Stufe tun, erhalten Sie ein völlig anderes Ergebnis (Kashima wird viel schlechter sein). Bis zum letzten ist es genau so, wenn es aus den Daten der 1. Stufe berechnet wird, also bin ich mir nicht sicher.

Wenn Sie darüber nachdenken, ist es natürlich, aber es liegt nahe an der 1. Etappe 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 Insbesondere scheint es eine allgemeine Korrelation zwischen "Gesamtpunkten und Gesamtzielen" und "Angriffs- und Verteidigungskraft" zu geben (obwohl es ein Problem wäre, wenn es keine gäbe).

schließlich

Ich habe versucht zu simulieren, was passieren würde, wenn ich den Toto der 2. Stufe mit diesem Vorhersagemodell kaufen würde, aber ich habe aufgehört, weil es ein schreckliches Ergebnis zu sein scheint (lacht). Es scheint, dass wir gute Vorhersagen treffen können, wenn die Wettbewerbsdichte in kurzer Zeit hoch ist. Sagen Sie 17 Verse mit Daten bis zu 16 Versen voraus. Ist das etwas besser?

Es ist erstaunlich, dass diese Art des Modellierens einfach geworden ist und viel Spaß macht.