I have a model that produces the expected points of every football (soccer) team for the following season, given the results from the previous season (the model considers the number of goals scored/conceded by each team when playing at home/away, using a Poisson distribution to determine the outcome of every game in the next season). I also have the data for the actual points each team achieved in the second season.

Given the predicted points for the next season, and the actual points, is it possible to calculate how likely each team is to finish in each position, out of 28 teams? (i.e. likelihood of coming 1st, 2nd, ..., 28th)

My initial thought was to create a range of points for each table position (i.e. 130+ points for first, 120-129 points for second, so on) and calculate the likelihood of a team falling within this range. I don't know how to implement this however.

Any replies would be much appreciated!

A portion of the table of actual vs predicted points is below:

enter image description here


1 Answer 1


I think this could be done using the parametric bootstrap.

Since you have a model for the expected number of points, and you make assumptions about the distribution of those points, you can essentially simulate new data. The approach is as follows (I'll demonstrate using some R code).

You have your data that you used to fit your model:

# Simulating some data to use.
nteams = 50
x = rnorm(nteams)
lambda = 100*exp(-0.25*x)
y = rpois(nteams, lambda)

fit=glm(y~x, family = poisson())

From your model, for each team, you can predict the expected number of points:

predicted_lambda = predict(fit, type='response')

Do the following 1000 times: For each team, simulate a poisson random variable using the expected points from your model.

nsims = 1000
# Store simulations here
sim_results = matrix(rep(0, nteams*nsims), nrow = nsims)

for(i in 1:nsims){
  #Simulate new data based on the fit
  sim_results[i, ] = rpois(nteams, predicted_lambda)


For each simulation you just did, order each team from first to last

for(i in 1:nsims){
   # This is equivalent to argsort in something like numpy
  sim_results[i, ] = sort(sim_results[i,], index.return=T)$ix

You now have 1000 simulations of the ordering of each teams. All you have to do now is count the rankings each column. For example, here is the probability that one of my simulated teams winds up in the indicated spot

sim_results[, 1] %>% 
  table %>% 

   3     6     7     8    12    14    19    20    25    28    35    39    43    50 
0.012 0.025 0.036 0.012 0.010 0.019 0.117 0.054 0.030 0.002 0.123 0.003 0.556 0.001 

In practice, you should also simulate the coefficients from their sampling distribution and use them to generate new predictions in each simulation step. I've spoken a little about this in this answer. The best approach, however, would be to use a Bayesian poisson regression (assuming the model is well specified) because you can integrate over all the uncertainty in the model.

  • 1
    $\begingroup$ @codemachino Can you share your code and a link to the data so I can check? $\endgroup$ Commented Oct 22, 2021 at 13:46
  • $\begingroup$ @codemachino Just add a link to your post. $\endgroup$ Commented Oct 22, 2021 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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