# Likelihood of a predicted value falling within a range of values

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: 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).

# 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 %>%
prop.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.

• @codemachino Can you share your code and a link to the data so I can check? Oct 22 at 13:46