Découvrez à quel âge un joueur de baseball peut frapper une base avec Python et Stan
Aperçu
Il y a beaucoup de choses à craindre dans le monde, mais je pense que la chose la plus importante est: "Quel âge ont les joueurs de baseball quand ils peuvent frapper le plus à domicile?" Cette fois, j'aimerais calculer le taux de réussite au domicile pour l'âge de chaque frappeur par une estimation bayésienne à partir des données des ligues majeures. La vraie définition est différente, mais par souci de simplicité, nous la définissons comme "taux de succès de base = hits / hits de base". En d'autres termes, le taux de réussite de la base est la valeur de la probabilité qu'un joueur frappe une base lorsqu'il se tient sur le siège (sauf dans le cas de quatre balles).
L'outil à utiliser est
- Python
- Pandas
- Stan (PyStan)
est.
Les lecteurs visés sont ceux qui s'intéressent ou étudient une partie du baseball, des statistiques bayésiennes, Python et Stan (presque personne, et vice versa, vous n'avez pas besoin de lire cet article. Il y a une inquiétude). Par conséquent, je ne vais pas expliquer chacun en détail, mais je voudrais énumérer quelques bonnes publications, le cas échéant.
Excuse
Je suis un profane de la statistique et un profane bayésien et un profane de Stan, donc le contenu suivant peut contenir divers types d'erreurs, de la plus fondamentale à la plus subtile. Je vous remercie(?).
Obtenez des données
La base de données de Lahman est un site incontournable pour tous les joueurs des ligues majeures au cours des 100 dernières années. , Le nombre de hits de la base d'accueil ...) est publié, je vais donc l'utiliser. Cette fois, nous allons télécharger la version 2016 du fichier csv ("Télécharger la version 2016" -> "2016 - version délimitée par des virgules").
Si vous décompressez le fichier téléchargé, vous perdrez la motivation car il y a des fichiers csv infinis, mais cette fois je vais l'utiliser
--Master.csv (informations de base des joueurs) --Batting.csv
Seulement.
Prétraitement des données
ʻImportla bibliothèque à utiliser cette fois, lisez le fichier csv et faites-lepandas.DataFrame`. Pandas est une bibliothèque de type Excel incroyablement pratique qui peut être utilisée sur Python (approprié), et est très utile pour conserver des données sous la forme d'un tableau et effectuer des opérations détaillées sur les données.
import numpy as np
import pandas as pd
import pystan
import matplotlib.pyplot as plt
bat = pd.read_csv('Batting.csv')
mas = pd.read_csv('Master.csv')
Une colonne dans Batting.csv est la performance d'un an avec un joueur. Voyons à quoi ça ressemble.
bat.tail()
Je vois, ça fait du bien. Différentes notes sont liées à l'année et à l'ID du joueur. Cependant, ce que je veux savoir cette fois, c'est la corrélation entre l'âge et les notes, je dois donc ajouter une colonne d'âge. En fait, dans Master.csv, l'identifiant du joueur et l'année de naissance sont alignés, je vais donc l'utiliser. Jetons un coup d'œil à la forme de Master.csv.
mas.tail()
D'ACCORD. Fusionnez les deux DataFrame en utilisant l'identifiant du joueur comme clé. Après cela, l'âge du joueur est calculé en soustrayant l'année de naissance du joueur de l'année de la saison. Cependant, dans les ligues majeures, l'âge des joueurs d'une saison donnée est officiellement utilisé à partir du 30 juin, donc pour les joueurs nés en juillet ou après, l'année de naissance devrait être l'année suivante.
mas['birthYear'] = mas['birthYear'] + (mas['birthMonth'] >= 7)
bat = pd.merge(bat, mas[['playerID', 'birthYear']], on='playerID')
bat['age'] = bat['yearID'] - bat['birthYear']
En passant, je voudrais mentionner les deux livres suivants qui contiennent des conseils pour faire des statistiques de baseball telles que l'histoire de l'âge au-dessus et divers résultats d'analyse.
Analysons-le, mais il semble que ce n'est pas si bon d'utiliser des données trop anciennes ou des données de joueurs qui n'ont participé que peu, donc cette fois je vais les couper de la base de données. Plus précisément, nous n'utiliserons que les données des joueurs nés après 1965 et disposant de plus de 1000 sièges au total. Pour ce faire, nous calculons d'abord le nombre de sièges (PA), puis utilisons groupby pour obtenir le nombre total de sièges pour ce joueur.
bat['PA'] = bat['AB'] + bat['BB'] + bat['HBP'] + bat['SH'] + bat['SF']
pa_c = bat.groupby('playerID').sum()[['PA']]
pa_c.columns = ['PA_Career']
bat = pd.merge(bat, pa_c, left_on='playerID', right_index=True)
Jetons un coup d'œil au DataFrame une fois.
bat.head()
Une colonne d'âge des joueurs a été ajoutée à droite, ce qui est bien. Maintenant, faisons une véritable analyse.
Modèle 0 (moyenne) et pourquoi une estimation bayésienne est nécessaire
Ci-dessous, nous allons commencer avec un modèle simple et le mettre à jour progressivement pour une meilleure sensation.
Tout d'abord, je voudrais essayer d'additionner le nombre de sièges et de bases de tous les joueurs pour chaque âge et de le diviser. Par exemple, "Si le nombre total de coups sûrs d'un joueur de 30 ans inclus dans les données est de 100 000 et que 3000 coups sûrs sont effectués à domicile, le taux de réussite à domicile d'un joueur de 30 ans est de 3 000/100 000, soit 0,03." Faire. Ce type de processus est très facile à faire avec la fonction «group by» de «Pandas».
HRdf = bat.groupby('age').sum()[['AB', 'HR']]
HRdf['HRR'] = HRdf['HR'] / HRdf['AB']
fig, ax = plt.subplots()
HRdf['HRR'].plot(kind='line', style='o-', c='black', ax=ax)
ax.set_ylabel('p')
L'axe horizontal représente l'âge et l'axe vertical le taux de réussite de base. Je vois. On peut voir que le taux de réussite à domicile augmente lentement à partir de 19 ans et est le plus élevé à 29 ans. La valeur à ce moment est de 0,034, soit environ 3,4%. Après cela, il a lentement diminué et après l'âge de 40 ans, il a fortement diminué.
Donc, il est normal de se retrouver avec l'idée que vous voulez connaître le taux de réussite de la base pour chaque âge, mais il semble y avoir un problème à considérer le graphique ci-dessus comme le «taux de réussite de la base réelle».
(a). ** Les différences individuelles ne sont pas prises en compte. ** Dans ce chiffre, les joueurs âgés de 37 et 8 ans sont censés frapper autant les bases que les joueurs à la fin de la vingtaine. Est-ce vrai? Si vous pensez calmement, il semble que ce soit parce que lorsque vous êtes dans la trentaine, les frappeurs avec des performances de frappe médiocres se retirent, et seuls les frappeurs qui peuvent frapper les bases restent. Ces frappeurs qui peuvent être actifs jusqu'à la fin de la trentaine auraient dû frapper plus de bases à la fin de la vingtaine, mais ce n'est pas clair sur ce chiffre. En effet, les différences individuelles ne sont pas prises en compte car tous les frappeurs ont été moyennés.
(b). ** Bruyant. ** Dans la figure ci-dessus, un frappeur de 39 ans frappe anormalement un port d'attache, mais cela ne signifie pas que la capacité d'un joueur de baseball à frapper un port d'attache augmente soudainement à l'âge de 39 ans. .. Je ne vois aucune raison pour laquelle je pourrais soudainement frapper un port d'attache simplement parce que j'avais 39 ans. En pensant normalement, les véritables capacités d'un joueur de baseball de 39 ans devraient être à peu près identiques ou légèrement inférieures à celles de 37 ou 38 ans.
(c). ** Je ne connais pas la fiabilité des points de données. ** Dans les données ci-dessus, les nombres de frappeurs de 30 ans et de frappeurs de 45 ans sont complètement différents. Dans le cas d'un traitement qui prend la moyenne globale comme cette fois, la fiabilité de la valeur augmente à mesure que le nombre d'échantillons augmente, de sorte que les données vieilles de 45 ans sont presque peu fiables (représentant une valeur proche de la valeur réelle). Il est peu probable que vous y soyez). Néanmoins, j'estime qu'il est étrange de s'aligner avec dignité sur les données de 30 ans.
Pour les raisons ci-dessus, j'ai envie de faire de la modélisation statistique. Ce qui précède est mon idée (amateur), donc si vous voulez apprendre les avantages et les idées de base de la modélisation statistique, par exemple, vous devriez lire le livre suivant.
- [Introduction à la modélisation statistique pour l'analyse des données - Modèle linéaire généralisé / Modèle hiérarchique de Bayes / MCMC](https://www.amazon.co.jp/%E3%83%87%E3%83%BC%E3] % 82% BF% E8% A7% A3% E6% 9E% 90% E3% 81% AE% E3% 81% 9F% E3% 82% 81% E3% 81% AE% E7% B5% B1% E8% A8 % 88% E3% 83% A2% E3% 83% 87% E3% 83% AA% E3% 83% B3% E3% 82% B0% E5% 85% A5% E9% 96% 80 __% E4% B8% 80 % E8% 88% AC% E5% 8C% 96% E7% B7% 9A% E5% BD% A2% E3% 83% A2% E3% 83% 87% E3% 83% AB% E3% 83% BB% E9 % 9A% 8E% E5% B1% A4% E3% 83% 99% E3% 82% A4% E3% 82% BA% E3% 83% A2% E3% 83% 87% E3% 83% AB% E3% 83 % BBMCMC-% E7% A2% BA% E7% 8E% 87% E3% 81% A8% E6% 83% 85% E5% A0% B1% E3% 81% AE% E7% A7% 91% E5% AD% A6-% E4% B9% 85% E4% BF% 9D-% E6% 8B% 93% E5% BC% A5 / dp / 400006973X / ref = pd_sim_14_9? _Encoding = UTF8 & psc = 1 & refRID = TCFDQTQGKQEM42AERGJF)
- Bayesian Methods for Hackers
Modèle 1 (Bayes hiérarchiques compte tenu de l'âge)
Faisons le modèle suivant.
J'ai soudainement émis une formule incompréhensible, je vais donc l'expliquer en japonais. Supposons que les coups de base d'un frappeur dans une saison suivent une distribution binomiale déterminée par les coups $ AB $ et le taux de succès de la base de base $ p $. De plus, supposons que le taux de réussite de la base de départ $ p $ soit représenté par une fonction logistique qui prend la somme de la constante $ \ beta $ et de la variable $ r_ {age} $ déterminée par l'âge du joueur comme argument. Ici, la valeur de $ r_ {age} $ pour chaque âge suit une distribution normale avec une moyenne de 0 et une variance de $ s_ {age} $. À partir de maintenant, nous utiliserons PyStan, qui est un wrapper pour Python du logiciel libre Stan que MCMC peut faire, pour effectuer une estimation bayésienne de chaque paramètre.
Ici, qu'est-ce que l'estimation bayésienne, qu'est-ce que MCMC, comment utiliser Stan, quelle est la distribution de probabilité, pourquoi les gens vivent en premier lieu, comment la conscience humaine est née, l'IA Je ne vais pas l'expliquer ici car il fera complètement noir et il y a beaucoup d'explications sur des gens plus décents que moi si je dis si je peux avoir conscience.
Pour les livres relativement faciles à lire en japonais, les bases de l'estimation bayésienne et de la modélisation statistique sont énumérées ci-dessus ["Introduction à la modélisation statistique pour l'analyse des données"](https://www.amazon. co.jp/%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E3%81%AE%E3%81%9F%E3%82 % 81% E3% 81% AE% E7% B5% B1% E8% A8% 88% E3% 83% A2% E3% 83% 87% E3% 83% AA% E3% 83% B3% E3% 82% B0 % E5% 85% A5% E9% 96% 80 __% E4% B8% 80% E8% 88% AC% E5% 8C% 96% E7% B7% 9A% E5% BD% A2% E3% 83% A2% E3 % 83% 87% E3% 83% AB% E3% 83% BB% E9% 9A% 8E% E5% B1% A4% E3% 83% 99% E3% 82% A4% E3% 82% BA% E3% 83 % A2% E3% 83% 87% E3% 83% AB% E3% 83% BBMCMC-% E7% A2% BA% E7% 8E% 87% E3% 81% A8% E6% 83% 85% E5% A0% B1% E3% 81% AE% E7% A7% 91% E5% AD% A6-% E4% B9% 85% E4% BF% 9D-% E6% 8B% 93% E5% BC% A5 / dp / 400006973X / ref = pd_sim_14_9? _encoding = UTF8 & psc = 1 & refRID = TCFDQTQGKQEM42AERGJF), mais pour Stan ["Modélisation statistique Bayes avec Stan et R"](https://www.amazon.co.jp/exec/obidos/ASIN/4320112423/statmodeling Si vous lisez -22 /), vous pouvez le comprendre.
Le code pour PyStan est ci-dessous.
AB = bat['AB'].values
HR = bat['HR'].values
age_map = {age: idx+1 for idx, age in enumerate(np.unique(bat['age'].values))}
age = bat['age'].map(age_map).values.astype(int)
data = dict(N=len(AB), N_age=len(np.unique(age)), HR=HR, AB=AB, age=age)
model_code = """
data {
int N;
int N_age;
int HR[N];
int AB[N];
int age[N];
}
parameters {
real beta;
real<lower=0> s_age;
vector[N_age] r_age;
}
transformed parameters {
vector[N] p;
for (i in 1:N)
p[i] = inv_logit(beta + r_age[age[i]]);
}
model {
r_age ~ normal(0, s_age);
HR ~ binomial(AB, p);
}
generated quantities {
real p_age[N_age];
for (i in 1:N_age)
p_age[i] = inv_logit(beta + r_age[i]);
}
"""
fit = pystan.stan(model_code=model_code, data=data, chains=4, iter=1000)
Lorsque vous exécutez le code ci-dessus, l'ordinateur fera de son mieux pour obtenir la distribution des différents paramètres ($ \ beta, s_ {age}, r_ {age} $). Dans l'estimation bayésienne, la valeur du paramètre n'est pas fixée à une valeur, mais est obtenue sous forme de distribution, il est donc possible d'exprimer l'incertitude des données. Par exemple, les données de 45 ans ont moins d'échantillons et sont plus incertaines que les données de 30 ans, ce qui se traduit par une distribution plus large des paramètres. En d'autres termes, le problème (c) de l'analyse simplement en faisant la moyenne des données expliquées ci-dessus sera résolu.
Après tout, je veux estimer le taux de réussite de la base à domicile $ p $ pour chaque âge, donc je vais le tracer avec $ r_ {age} $ avec le code suivant.
def plot_r_p(fit):
inv_age_map = {v: k for k, v in age_map.items()}
x_age = [inv_age_map[a] for a in np.unique(age)]
r_age_trace = fit.extract()['r_age']
p_age_trace = fit.extract()['p_age']
#Médiane et 50%・ 95%Calculez le centile aux points suivants pour tracer l'intervalle prédit
percents = [2.5, 25, 50, 75, 97.5]
r_age = []
p_age = []
for percent in percents:
r_age.append(np.percentile(r_age_trace, percent, axis=0))
p_age.append(np.percentile(p_age_trace, percent, axis=0))
fig, axs = plt.subplots(nrows=1, ncols=2, sharex=True, figsize=(12, 4))
# axs[0] : r_age
axs[0].plot(x_age, r_age[2], 'o-', color='blue') #Médian
axs[0].fill_between(x_age, r_age[0], r_age[4], alpha=0.1, color='blue') # 95%section
axs[0].fill_between(x_age, r_age[1], r_age[3], alpha=0.3, color='blue') # 50%section
axs[0].set(xlabel='age', ylabel=r'$r_{age}$')
# axs[1] : p_age
axs[1].plot(x_age, p_age[2], 'o-', color='red')
axs[1].fill_between(x_age, p_age[0], p_age[4], alpha=0.1, color='red')
axs[1].fill_between(x_age, p_age[1], p_age[3], alpha=0.3, color='red')
axs[1].set(xlabel='age', ylabel=r'$p_{age}$')
return fig, axs
fig, axs = plot_r_p(fit)
HRdf['HRR'].plot(kind='line', style='o-', c='black', ax=axs[1]
Les données bleues sur la gauche sont $ r_ {age} $, les données noires sur la droite sont la moyenne simple indiquée ci-dessus et les données rouges sont le taux de réussite estimé à domicile pour chaque âge $ p_ {age} $. Dans les deux graphiques, la couleur de fond clair représente la section prédite à 50% et la couleur d'arrière-plan sombre représente la section prédite à 95%.
Donc, quand vous regardez les résultats, vous avez l'impression que hmm. Je pense que le problème (c) a été partiellement résolu, mais rien n'a été fait pour (a) les différences individuelles et (b) les problèmes de bruit. Par exemple, comme d'habitude, à l'âge de 39 ans, il est censé frapper soudainement une base d'attache, et comme nouveau problème, sa capacité s'améliore après la quarantaine, qui a peu de points de données, ce qui est un résultat irréaliste (ceci). Je pense que c'est probablement parce que la probabilité est plus élevée lorsqu'elle se rapproche de la moyenne de tous les âges).
Modèle 2 (compte tenu des séries chronologiques)
Afin de résoudre les problèmes ci-dessus, nous examinerons la série chronologique. Cela signifie que le modèle précédent traitait les capacités spécifiques à l'âge de manière complètement indépendante (correspondant à la distribution normale de $ r_ {age} $ avec une moyenne de 0 dans l'équation du modèle 1). Cependant, si vous y réfléchissez, votre capacité à l'âge j devrait être proche de votre capacité à l'âge (j-1). Si cela se reflète dans le modèle, je pense que ce sera comme suit.
Dans ce modèle, la moyenne de la distribution normale que $ r_ {age} [j] $ suit est $ r_ {age} [j-1] $. Cela représente l'hypothèse que "la capacité à l'âge j devrait être proche de la capacité à l'âge (j-1)".
AB = bat['AB'].values
HR = bat['HR'].values
age_map = {age: idx+1 for idx, age in enumerate(np.unique(bat['age'].values))}
age = bat['age'].map(age_map).values.astype(int)
data = dict(N=len(AB), N_age=len(np.unique(age)), HR=HR, AB=AB, age=age)
model_code = """
data {
int N;
int N_age;
int HR[N];
int AB[N];
int age[N];
}
parameters {
real beta;
real<lower=0> s_age;
vector[N_age] r_age;
}
transformed parameters {
vector[N] p;
for (i in 1:N)
p[i] = inv_logit(beta + r_age[age[i]]);
}
model {
r_age[1] ~ normal(-sum(r_age[2:N_age]), 0.001);
r_age[2:N_age] ~ normal(r_age[1:(N_age-1)], s_age);
HR ~ binomial(AB, p);
}
generated quantities {
real p_age[N_age];
for (i in 1:N_age)
p_age[i] = inv_logit(beta + r_age[i]);
}
"""
fit2 = pystan.stan(model_code=model_code, data=data, chains=4, iter=1000)
Le seul changement par rapport au modèle 1 est la partie où r_age [2: N_age] ~ normal (r_age [1: (N_age-1)], s_age) est défini dans model (somme de $ r_ {age} $). Je pense que la partie où r_age [1] ~ normal (-sum (r_age [2: N_age]), 0.001) est mis à 0 n'est vraiment pas très bonne ...).
Le résultat est tracé comme suit.
Je vois ~. Il y a deux avantages par rapport au modèle 1.
――L'augmentation des capacités à l'âge de 39 ans a été un peu supprimée. ――La capacité dans ma quarantaine diminue d'année en année
Tous ces effets tiennent compte des séries chronologiques et sont considérés comme des résultats réalistes. Donc, je pense que le bruit du problème (b) s'est un peu amélioré. Alors, enfin, faisons un modèle qui tient compte des différences individuelles.
Modèle 3 (compte tenu des différences individuelles)
Alors, ajoutez un terme $ r_ {id} $ qui représente les différences individuelles du modèle. C'est une variable qui représente la différence de capacité de chaque joueur, et devrait suivre une distribution normale avec une moyenne de 0.
id_map = {id: idx+1 for idx, id in enumerate(np.unique(bat['playerID'].values))}
id = bat['playerID'].map(id_map).values.astype(int)
data = dict(N=len(AB), N_age=len(np.unique(age)), N_id=len(np.unique(id)), HR=HR, AB=AB, id=id, age=age)
model_code = """
data {
int N;
int N_id;
int N_age;
int HR[N];
int AB[N];
int id[N];
int age[N];
}
parameters {
real beta;
real<lower=0> s_age;
real<lower=0> s_id;
vector[N_age] r_age;
vector[N_id] r_id;
}
transformed parameters {
vector[N] p;
for (i in 1:N)
p[i] = inv_logit(beta + r_age[age[i]] + r_id[id[i]]);
}
model {
r_id ~ normal(0, s_id);
r_age[1] ~ normal(-sum(r_age[2:N_age]), 0.001);
r_age[2:N_age] ~ normal(r_age[1:(N_age-1)], s_age);
HR ~ binomial(AB, p);
}
generated quantities {
real p_age[N_age];
for (i in 1:N_age)
p_age[i] = inv_logit(beta + r_age[i]);
}
"""
fit3 = pystan.stan(model_code=model_code, data=data, chains=4, iter=1000)
Voyons le résultat.
Je vois ~. N'est-ce pas un très bon résultat? L'âge change en douceur avec chaque capacité et semble assez raisonnable. Ici, la valeur estimée de $ p_ {age} $ est beaucoup plus petite que la valeur mesurée, car la valeur estimée représente "$ p_ {age} $ du joueur moyen par âge". Je pense. Je pense que plus les joueurs sont incroyables, plus ils auront de sièges dans la mesure réelle, donc si vous prenez la moyenne globale, le poids des joueurs incroyables sera plus élevé que celui des joueurs qui ne sont pas très lourds, et par conséquent, la valeur sera supérieure au $ p_ {age} $ du joueur moyen. Je vais sortir.
En conclusion, les joueurs de 29 ans sont les plus susceptibles de frapper les bases, un peu plus de 2,5% pour le joueur moyen. En regardant les changements de capacités avant et après l'âge de 29 ans, vous pouvez voir que la croissance avant l'âge de 29 ans est plus rapide que la baisse après l'âge de 29 ans. De plus, par exemple, $ p_ {age} $ autour de 20 ans est un peu plus élevé que 1,5%, donc en moyenne, les joueurs d'environ 30 ans ont $ (2,5 / 1,5) \ simeq 1,7 que les joueurs d'environ 20 ans. Vous pouvez vous attendre à toucher une base environ $ fois.
Je pense que je peux encore faire diverses choses, mais pour le moment, par ici ...
Résumé
Ainsi, à la suite de l'estimation bayésienne utilisant PyStan, l'âge auquel un joueur de baseball peut frapper le plus souvent la base, il a été constaté que le joueur de 29 ans peut le plus frapper la base.
c'est tout, merci beaucoup.
Bonus: j'ai recherché d'autres indicateurs de la même manière, mais je pense que les énumérer est hors de portée de qiita, donc c'est différent (http://hoture6.hatenablog.com/entry/2017/06). / 08/212256). Si vous êtes intéressé, s'il vous plaît.
Recommended Posts