(Bundesliga) game prediction¶
Simple model for the predictions of soccer games: How many goals a team scores.
As data only the results from prior games are used.
In [6]:
#top_games
Poisson distribution¶
Probability for outcome $k \in \{0, 1, 2, \dots\}$
$$ P_\lambda (k) = \frac{\lambda^k}{k!}\, \mathrm{e}^{-\lambda} $$with parameter $\lambda>0$
- $\lambda$ is also the expectation and variance of the distribution
In [20]:
import scipy.stats
k=np.arange(0,10)
lambda_= 3.2
plt.figure(figsize=(8,6))
plt.plot(k, scipy.stats.poisson.pmf(k, lambda_), 'bo', ms=6, label='poisson pmf')
plt.xlabel("k")
plt.ylabel("probability mass")
scipy.stats.poisson.pmf(k, lambda_)
Out[20]:
Probabilistic Model¶
Each team $i$ has a offence and defence strength (distribution). (Note that the average goals per game $\approx 3 \Rightarrow \Delta \mu=1.5$):
$$ offence_i \sim \mathcal N(\mu=1.5, \tau=1) $$$$ defence_i \sim \mathcal N(\mu=0, \tau=1) $$$\mathcal N$ is the Gaussian distribution with parameters
- mean: $\mu$
- precision: $\tau=1/\sigma^2$ (variance: $\sigma^2$)
Model: The number of goals that team $i$ scores against team $j$ is Poisson distributed with
$$ goals_{ij} = Poisson \left(\lambda = (offence_i-defence_j) \right) $$Graphical representation of the model¶
In [24]:
plot_model()
Implementation with pymc¶
In [43]:
offence = [pymc.Normal("offence%i"%i, tau=1., mu=1.5) for i in top_club_ids]
defence = [pymc.Normal("defence%i"%i, tau=1., mu=0.) for i in top_club_ids]
home_goals = [pymc.Poisson("home_goals_%i_%i" % (heim,gast),
observed=True, value=goals,
mu = pymc.Lambda("home_rate_%i_%i" % (heim, gast),
lambda o=offence[heim], d=defence[gast]: low if o-d<0. else o-d))
for i,(heim, gast, goals) in enumerate(home_goals_)]
away_goals = [pymc.Poisson("away_goals_%i_%i" % (heim,gast),
observed=True, value=goals,
mu = pymc.Lambda("away_rate_%i_%i" % (heim, gast),
lambda o=offence[gast], d=defence[heim]: low if o-d<0. else o-d))
for i,(heim, gast, goals) in enumerate(away_goals_)]
Sampling with pymc¶
In [45]:
mcmc = pymc.MCMC([offence, defence, home_goals, away_goals])
nb_samples=12000
burn=4000
mcmc.sample(nb_samples, burn=burn)#, thin=50)
In [52]:
# sampling of probable results according to the model
#print_sample_results(most_probable=False)
# Sample results from a sample of the trace
nb_p, p, _ = print_sample_results(most_probable=True)
Sampling histograms¶
In [56]:
nb_clubs = top_club_ids.max() + 1
bins=40
fig, axes = plt.subplots(nrows=nb_clubs, ncols=2, figsize=(10, 50))
for i in top_club_ids:
title = "Offence of " + clubs[clubs.index==i]["Name"][i].decode('latin')
axes[i, 0].set_title(title)
axes[i, 0].hist(mcmc.trace("offence%i"%i)[:], bins=bins, range=(0,4.2))
axes[i, 1].hist(mcmc.trace("defence%i"%i)[:], bins=bins, range=(-2.,2.2))
title = "Defence of " + clubs[clubs.index==i]["Name"][i].decode('latin')
axes[i, 1].set_title(title)
#fig.suptitle("Offence and defence distribution of the clubs.")
fig.subplots_adjust(hspace=0.5)
fig.tight_layout()
Exercise: Distribution of expected goals¶
Use the model and the mcmc trace to predict how many goals a teams scores agains anouther team. Programm the corresponding python (plot) functions, e.g.
In [64]:
# id 0: Bayern München
# id 17: SV Darmstadt
plot_goal_diffs(0, 17)
# probability that Bayern München scores 0,1,2, ... goals against SV Darmstadt
# array([ 6.71642520e-02, 1.75088332e-01, 2.34376889e-01, ...
(Possible) Solution¶
In [32]:
# Hierachical model
offence_mean = pymc.Normal("offence_mean", mu=1.5, tau=1.)
offence = [pymc.Normal("offence%i"%i, mu=offence_mean, tau=1.) for i in top_club_ids]
defence = [pymc.Normal("defence%i"%i, mu=0., tau=1.) for i in top_club_ids]
#midfield = [pymc.Normal("midfield%i"%i, mu=1., tau=1.) for i in top_club_ids]
om = offence #[offence[i] * midfield[i] for i in top_club_ids]
dm = defence #[defence[i] * midfield[i] for i in top_club_ids]
home_advantage = ha = pymc.Normal("home_advantage", tau=10., mu=1.0)
#noise_prior_g = pymc.HalfNormal("noise_prior_g", tau=1.)
#noise_h = [pymc.Normal("noise_h_%i_%i"% (heim, gast), mu=0, tau=noise_prior_g)
# for (heim, gast, goals) in home_goals_]
#noise_a = [pymc.Normal("noise_a_%i_%i"% (heim, gast), mu=0, tau=noise_prior_g)
# for (heim, gast, goals) in away_goals_]
home_goals = [pymc.Poisson("home_goals_%i_%i" % (heim,gast),
observed=True, value=goals,
mu = pymc.Lambda("home_rate_%i_%i" % (heim, gast),
lambda o=om[heim]*ha, d=dm[gast], n=0.: low if o-d+n<0. else o-d+n))
for i, (heim, gast, goals) in enumerate(home_goals_)]
away_goals = [pymc.Poisson("away_goals_%i_%i" % (heim,gast),
observed=True, value=goals,
mu = pymc.Lambda("away_rate_%i_%i" % (heim, gast),
lambda o=om[gast], d=dm[heim]*ha, n=0.: low if o-d+n<0. else o-d+n))
for i, (heim, gast, goals) in enumerate(away_goals_)]
### total result (tie, home win, home loose)
noise_prior = pymc.HalfNormal("noise_prior", tau=1.)
noise = [pymc.Normal("noise_%i_%i"%(heim,gast), tau=1./(noise_prior**2))
for i, (heim, gast, goals) in enumerate(home_goals_)]
defence_scale = pymc.HalfNormal("defence_scale", tau=1.)
skills = [(offence[i] + defence_scale * defence[i]) for i in top_club_ids]
skills_difference = [pymc.Lambda("skill_difference_%i_%i" % (heim, gast),
lambda ho_s=skills[heim]*ha, aw_s=skills[gast], n=noise[i]:ho_s-aw_s+n)
for i, (heim, gast, result) in enumerate(results_)]
params = pymc.HalfNormal("params", tau=1., size=(3))
pymc.abs(skills_difference[i])
# soft_max outputs
draw_t = [pymc.exp(params[0] - params[1] * pymc.abs(skills_difference[i]))
for i in range(len((results_)))]
home_win_t = [pymc.exp(skills_difference[i] * params[2])
for i in range(len((results_)))]
away_win_t = [pymc.exp(- skills_difference[i] * params[2])
for i in range(len((results_)))]
sm_denominator = [draw_t[i]+home_win_t[i]+away_win_t[i] for i in range(len((results_)))]
p_ = [(draw_t[i]/sm_denominator[i], home_win_t[i]/sm_denominator[i], away_win_t[i]/sm_denominator[i])
for i in range(len((results_)))]
prediction = [pymc.Categorical("prediction_%i_%i" % (heim, gast),
observed=True, value=result,
p=p_[i])
for i, (heim, gast, result) in enumerate(toto_results_)]
#prediction_ = [pymc.Categorical("prediction_%i_%i_" % (heim, gast),
# p=p_[i])
# for i, (heim, gast, result) in enumerate(toto_results_)]
In [34]:
mcmc = pymc.MCMC([offence_mean, offence, defence, home_goals, away_goals, home_advantage, params, prediction, noise_prior, noise, defence_scale])
nb_samples=10000
burn=2000
mcmc.sample(nb_samples, burn=burn)
In [35]:
mcmc.trace("defence_scale")[:].mean(axis=0)
Out[35]:
Offence¶
The offence (mean) shift corresponds to the number of goals per game (times 2):
In [37]:
plt.hist(mcmc.trace("offence_mean")[:], bins=20)
plt.title("Offence mean distribution")
Out[37]:
Home advantage¶
In [38]:
plt.hist(mcmc.trace("home_advantage")[:], bins=20)
plt.title("Home advantage distribution")
Out[38]:
In [39]:
#mcmc.trace("defence0")
#offence