0
$\begingroup$

I have a general question about Bayesian inference which may help me solve a problem I have. It is best to illustrate this with an example. Inspired from this great post by AllenDowney:

https://github.com/AllenDowney/BiteSizeBayes/blob/master/08_soccer_soln.ipynb

Let's say I want to get an estimate for the number of goals scored by a team in a football match. I model this as a Poisson process with parameter $\lambda$. On $\lambda$ I put a Gamma distribution as the prior, with parameter $\alpha$. Let's say I use a prior value of 1.4 for $\alpha$ (avg number of goals scored by a team).

In the blog post, we are computing the posterior after a match has been played. So, given a game where 4 goals were scored by a team, we compute the posterior for that team which is now shifted to the right.

What would we do if we wanted to update the estimate in real-time? So instead of computing the posterior after the game has been played, we do this every 5 minutes until we hit 90 minutes. So I am interested in getting a posterior for 90 minutes given the data after 5 minutes, 10 minutes, 15 minutes, etc.. I can think of two ways of doing this. Please help me understand what method, if any, makes sense:

  1. We propagate the data we have to the expected number of goals in 90 minutes after $x$ minutes have passed. So, if after 30 minutes a team scored 1 time, we propagate that to (90/30)*1=3 goals and use this as input to compute the posterior using the same process as in the blog post. Issue: I feel like this is not correct, because what if a team scores in the first 5 minutes? Doing this would mean we expect the team to score 18 (!) times in a match. The first posterior calculation after 5 minutes would be complete trash. Although this should eventually converge to a realistic value as we do more updates?

  2. We don't use a value of 1.4 for $\alpha$ because that's based on a 90 minute match. Instead, we use alpha = 1.4 / (90/5) = 0.08 to scale it to a 5-minute prior. So, on average, a team scores 0.08 times every 5 minutes. We now do the posterior calculation every 5 minutes instead of every 90 minutes. Issue: I don't understand how we now get a prediction for the 90 minute posterior because the posterior we calculate every time will be based on 5 minutes. Also, how do we link the second posterior (for minutes 5-10) to the first (minutes 0-5)?

Perhaps I am missing something basic here. I really want to understand Bayesian inference better but feel like I am not completely getting it. Thanks!

$\endgroup$
2
  • $\begingroup$ does my post answer your question? $\endgroup$ – Demetri Pananos Apr 22 '20 at 4:47
  • $\begingroup$ Yes it's super helpful! Sorry for not getting back earlier. $\endgroup$ – BarkingCat Apr 23 '20 at 11:02
2
$\begingroup$

Ok, so from what I understand from the blog post...

  • The likelihood for the number of goals scored is Poisson. Each team has a goal scoring rate, $\lambda$ measured in units per game. We can divide this by 18 to yield the goal scoring rate per 5 minute increments.

  • A gamma prior is put on $\lambda$. This makes things particularly nice because the gamma prior is conjugate for the Poisson likelihood, so our posterior will be gamma as well. For some strange reason, $\beta=1$ in scipy's parameterization by default (see the documentation for the Gamma density here and contrast it with the Gamma density here. For consistency with the blog post, I will adopt the assumption that $\beta=1$).

I think the problem you are experiencing is that your likelihood is on the scale of 5 minutes, but the prior is on the scale of games. Not to worry, this should be an easy fix. Thanks to conjugacy, the posterior is gamma distributed as well. In particular, on the scale of games

$$ \lambda | y \sim \operatorname{Gamma}(1.4+ y_i, 1+1) $$

Here, $y_i$ is the number of goals scored in a single game. To change this to the scale of 5 minutes, we need to do some arithmetic.

Let $\tilde{n}$ be the number of 5 minute increments which have been observed. Then, a single game has 18 5 minute increments. So, our posterior should then look like

$$ \lambda | y \sim \operatorname{Gamma}(1.4+ \sum_i \tilde{y_i}, 1+ \dfrac{\tilde{n}}{18}) $$

Now, $\tilde{y}_i$ is the number of goals scored in the $i^{th}$ 5 minute increment. Let's make sure this is indeed out posterior. Here is some python code to double check. I use pymc3 to sample from the posterior and then compare against my analytical result.

import numpy as np
from scipy.stats import gamma, poisson
import matplotlib.pyplot as plt
import pymc3 as pm

x = np.linspace(0,10)
alpha = 1.4
# draw a lambda from the prior
lams = gamma(a = alpha).rvs(1)
# use this to draw 18 observations.  Each element is the number of goals scored
# by this team in the ith 5 minute increment
goals = poisson(mu=lams/18).rvs(size=18)

with pm.Model() as model:

    #This model is on the scale of games.
    lam = pm.Gamma('lam',alpha=1.4, beta = 1)
    y = pm.Poisson('y', mu=lam, observed=goals.sum())

    trace = pm.sample()

# Histogram of the posterior
plt.hist(trace['lam'], density = True)
# Plot the true parameter value in red
plt.axvline(lams, color = 'red')

# Now, plot the posterior using the analytical result on the scale of games.
# At the end of the game, this posterior should look like the one from pymc3
y = pm.Gamma.dist(alpha = 1.4 + sum(goals), beta = 1+1).logp(x).eval()

plt.plot(x, np.exp(y))

enter image description here

The posterior I computed analytically on the scale of 5 mins (orange) is very similar to posterior I obtained from pymc3 (histogram in blue) on the scale of games, which leads me to believe I am correct.

Here is a more concrete example. Suppose we observe the following over an entire game.

0, 0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 1, 2, 0, 0, 0, 0, 0

In the first 5 minutes, there are no goals. In the second 5 minutes, there are no goals, etc. How would our posterior look after each 5 minutes? Like so

enter image description here

Here is the code to perform the updating after each 5 minutes

fig, ax = plt.subplots(dpi = 120, ncols=6, nrows=3, figsize = (20,10), sharex = True, sharey = True)
ax = ax.ravel()

alpha = 1.4
beta = 1
num_goals = 0
num_5_mins = 0

for i, g in enumerate(goals):

    num_goals+=g
    num_5_mins +=1


    y = gamma(a = 1.4 + num_goals, scale = 1/(1+num_5_mins/18)).pdf(x)

    ax[i].plot(x,np.exp(y))
    ax[i].set_title(f'Goals Scored So Far: {num_goals}')
    ax[i].set_ylabel('Density')
    ax[i].set_xlabel(r'$\lambda$')


$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.