Performance prediction with a hierarchical poisson model using Template Model Builder

Eivind Baltzersen

• Master's thesis in Applied Physics and Mathematics
• Submission date: 2020 July
• Supervisor: Jarle Tufto

Norwegian University of Science and Technology
Department of Mathematical Sciences

Notation and terms

1 Introduction

Ranking is a way to compare different objects given their properties. It is ideally transitive, meaning that for $n$ objects, we may label them with $\mathit{i}\in [1..n]$ to order them. An important aspect of competitive sports and games is to rank teams and players in order to sufficiently determine who is to be labeled $1$, (also known as the winner).

In a sports context, the ranking often depends on a score, which is determined based on the performance of the team or player against another team or player. This is achieved through each team or player playing against other teams and players, for instance through a knockout or a knockout tournament. Examples from football¹ in Norway include the knockout cup NM i fotball, and the top level round robin tournament is called Eliteserien.

To determine the winner of a football match, we simply look at the team with the most goals within two 45 minute halves of a 90 minute match. Some tournaments also have overtime if the teams are tied after 90 minutes, but we will mainly study Eliteserien, which doesn't include overtimes. The winner of a match gets three points, the loser none. A tie awards both teams with one point. The winner of the league is the team with the most points.

We make note of earlier work, such as 1, which attempts an attack-defence model that we will propose here, but with a different hierarchical model. A simple bivariate poisson model (non-hierarchical) has also been attempted in 2. A model using the skellam distribution has also been attempted in 3. A comparison of a poisson scoring models with a goal shots model is compared in 4. There are many other papers for similar models, both in football and other sports and games competitions.

The ideas in this paper are mostly based on earlier works. The novel methods is applying TMB as a tool to implement the models, and comparing different time series models.

As a small disclaimer, I will mention that I have little understanding of football, and have not watched a single football match during the writing of this thesis. So any statements about football may be incorrect.

1.1 Tools

The main programming languages have been R and C++, with the TMB library. Results were stored in json format. Plots have been made using Python, with matplotlib, NumPy and SciPy. For more details

2 Theory

We will go through several statistical concepts; it will be assumed that the reader has some knowledge of linear algebra. We will not go through each theorem in depth, nor the derivation of them. For a comprehensive explanation, a textbook should be consulted.

• Distributions
• Models
• Ranking
• Template model builder (TMB)
• Prediction
• Quality measures

2.1 Distributions

2.1.1 Normal and log-normal distribution

The normal² distribution is a normal distribution, and is described by the density function in equation (1)

figure 1: Example normal distribution

The normal distribution is fully parameterised by its mean and variance:

With the (unbiased) estimators

where $c(n)$ is a bias-correction factor, because $E(s)<\sigma$. For large samples this is close to one, so this factor is typically ignored, but exact and approximate values for the normal distribution can be found. 5 6 7

The log-normal distribution is similar, with a change of variables $x\to \log (x)$:

We note that $log(Y)$ is normal distributed, which we will make use of later, as the normal distribution is simpler.

The multivariate form of the normal distribution, often abbreviated to MVN, is

with the parameters

$\Sigma$ is usually referred to as a covariance matrix.

We will consider a special parametrisation of the covariance matrix. The motivation is to removed restrictions on the entries. First of all, the covariance matrix is symmetric, so almost half of the entries are redundant when specifying the matrix. A more complex restriction, is that it must also be positive definite, i.e. ${(x-\mu )}^{\top }\Sigma (x-\mu )\ge 0$ and its diagonal entries are strictly positive.

First we consider the relation between the covariance matrix and the correlation matrix:

where ${\rm P}$ is the correlation matrix, which have the nice property that its diagonal consists of units. We can find a lower-triangular square root $L$, such that $P=L{L}^{\top }$. However, $L$ does not have a unit diagonal, for this we multiply by the inverse of its diagonal to obtain $\Theta =\overline{diag(L)}L$, as shown in equation (12)

We let $\mathbf{\theta }$ be the vectorisation of the lower triangular matrix, i.e. $\mathbf{\theta }=({\theta }_1,{\theta }_2,\dots ,{\theta }_{\mathit{k}+n})$.

We may also consider the elementwise root-log transform of the diagonal of $\Sigma$:

Thus, we may parametricise $\Sigma$ as from the vector $\mathbf{\theta }\oplus log(\mathbf{\sigma })\in {ℝ}^{T(\mathit{k})}$, where $T(\mathit{k})$ is the $\mathit{k}$-th triangular number by reversing the above steps. 8

2.1.2 Binomial distribution and the fundamental theorem of statistics

In determining the outcome of a match, we are either right or wrong. If we have a rule determining the correct outcome of a match with probability $p$, we have a binary distribution³.

Repeated determination of $n$ matches, results in a binomial distribution, where we expect the number of the correct guessed outcomes $m$, such that $m/n\approx p$.

The limiting distribution of a binomial variable $X_n$ as $n\to \infty$, can be approximated by a normal distribution. This theorem is sometimes referred to as the de moivre–laplace theorem, which is a special case of the central limit theorem⁴.

Typically, this approximation is practical for $n$ greater than $30$.

The estimator for $p$, $\widehat{p}$, is simply the number of correct guessed $m$ over the total outcomes. $m/n=\widehat{p}$. So the confidence interval, assuming normality, is of the form

We reiterate that this is unreliable for a small sample size, and other confidence intervals also exist. A list of other methods can be found in 9.

2.1.3 Poisson distribution

The poisson distribution is defined as the number of events occuring a fixed interval. It's described by the mass function in equation (17)

figure 2: Example poisson distribution

A remarkable property of the poisson distribution, is the simple relation between the mean and variance:

2.1.4 Skellam distribution

The difference between two independent poisson distributed variables is distributed by the skellam distribution. 10 The skellam distribution is of the form show in equation (19)

where $I_{\mathit{k}}$ is the modified bessel function of the first kind. The formula is complicated, but is shown below 11, page 375

2.1.5 Discrete $\mathrm{VAR}(1)$ model

The $\mathrm{VAR}(1)$ series is the multivariate generalization of the one-dimensional $\mathrm{AR}(1)$ series, which is a special case of the $\mathrm{AR}(p)$ series:

where ${\epsilon }_t\sim \mathrm{WN}(0,{\Sigma }_{\epsilon })$. 12, page 84 In our case, we consider $p=1$. We will also be assuming that $c=0$. The distribution is stable as long as $\left|{\varphi }_1\right|<1$.

figure 3: Example discrete VAR(1) model

The unconditional expectation and covariance are given by 13

and the unconditional distribution of $x_t$ is

The conditional expectation and covariance are given by

and the conditional distribution of $x_t$ on $x_{t-1}$ is given by

The multivariate version of the unconditional and conditional distribution is: 14

We will later be using equation (28) and equation (29), with zero-shifted mean.

2.1.6 Continuous $\mathrm{VAR}$ model

The continuous version⁵ is often attributed to Ornstein, Uhlenbeck and Vašíček. 15 16 The differential form and the formal solutions are shown in equation (30-31).

where the absolute value of the eigenvalues of $\theta$ should be strictly positive in the real part to be stable. 17, page 11 The relation between these two form can be found in the appendix, in equation (120-130)

figure 4: Example continuous VAR model

The unconditional expectation and covariance are given by 18

We will be assuming $\mu =0$, but we state the general forms for completeness.

and the unconditional distribution of $x_t$ is 19

The conditional expectation and covariance are given by 20

and the conditional distribution of $x_t$ on $x_{t-1}$ is given by 21, page 11

The multivariate version of the unconditional and conditional distribution is: 22

As with the $\mathrm{VAR}(1)$ process, we will later be using equation (38) and equation (39), with zero-shifted mean.

The relation to the $\mathrm{VAR}(1)$ model is not immediately obvious, but equality can be shown with the following substitutions: 2324, page 8

so the following two equations are equivalent

2.2 Change of variables in a density function

In some cases it may be useful to transform the variables, such as $Y=g(X)$, because we may have more tools available for the transformed density. Such as log-transforming a log-normal variable, to obtain a normal variable.

To be precise, this is only valid as long as $g$ is a strictly increasing function. For the decreasing case we add a sign to either side; the two cases can be unified by applying the absolute value to both sides. 25

and the relation between $f_Y$ and $f_X$ is given by

If we let $Y=exp(X)\sim \mathrm{Lognormal}(\mu ,{\sigma }^2)$, then

And for the opposite case, $X=log(Y)\sim \mathcal{N}(\mu ,{\sigma }^2)$, then

2.3 Likelihood function

The most general form of the likelihood function is defined as a mass/density function of a parameter $\theta$ given an outcome $x$

In many cases, $x$ is a tuple of i.i.d. variables, and $p_{\theta }(x)$ is a joint probability distribution of independent variables, and can be written as a product. So if $x$ is a vector of size $n$, we have:

2.3.1 Maximization and log-likelihood

The likelihood function represents the probability of obtaining $x$ for a given $\theta$. A reasonable assumption is then that $x$ is realized from a distribution where it has a high likelihood to be observed. So the $\theta$ that yields the highest likelihood is a natural candidate for determining the distribution of $x$. We seek to determine

Product are often cumbersome, and to simplify the above computation, we can apply the logarithm to $L$ to obtain a sum. Because the logarithm is a strictly increasing function, the maximum of $L$ is also the maximum of $log(L)$. We denote this logarithm by $\ell$:

From which we conclude that $x$ most likely derived from the distribution $p(x|\widehat{\theta })$.

2.4 Hierarchical modelling and empirical bayes method

Bayesian⁶ hierarchical modelling is a way to describe a hierarchy of distributions, where the parameters of the upper layers are dependent on the distributions of lower layers.

The simplest such model is the two-stage hierarchical model, as shown in equation (56-58)

where $Y$ is the observed data, $\theta$ is a parameter, and $\varphi$ is a hyperparameter. $p(\theta ,\varphi )$ is the prior distribution. In bayesian hierarchical modelling the hyperparameter is given a hyperprior, a distribution on the parameter, or $p(\varphi )$.

In Empirical Bayes the hyperparameter is a fixed value. The estimation of this hyperparameter will be found using MLE.

The posterior theorem or bayes theorem is shown below in equation (59)

The distribution for $Y$ can then be found by marginalization⁷ over the parameters, or random effects.

2.5 Bradley–terry model

The bradley–terry model 26 is used to make paired comparisons of individuals in a transitive way, using a single parameter.

Where $p_{\mathit{i}}=e^{{\beta }_{\mathit{i}}}$ and $\sigma (x)=\overline{1-exp(x)}$ is the logistic function (a sigmoid function). The relation between the $\mathrm{logit}(x)=log(\frac{x}{1-x})$ logit function with the logistic function is that they are inverses, i.e. $\mathrm{logit}(x)=\bar{\sigma }(x)$, so we also have the identity:

An example where this model is used, is in elo ranking, most commonly known from chess. A player's rating is given by a number $r_{\mathit{i}}$, which is related to ${\beta }_{\mathit{j}}$ by ${\beta }_{\mathit{i}}=log(10)/400·r_{\mathit{i}}$. So a player with $r_1=1700$ playing against a player with $r_2=1900$ will yield the following probabilities:

A win or a loss updates the elo rating of each player, but the bradley–terry model has no mechanism for updating. Nor does it tell us how to initially calculate ratings, which must be found with inference.

We also note that the bradley–terry model does not handle ties, only binary win/loss outcomes.

The winner is the team with the highest score in a match. If the score is poisson distribution, we can determine the result by checking the value of $y_{\mathit{i}\mathit{j}}^{\mathrm{Home}}-y_{\mathit{i}\mathit{j}}^{\mathrm{Away}}$, and checking if it is strictly positive, zero or strictly negative.

2.5.1 Skellam distribution and the bradley–terry model

By defining the best team as having the score $p_M=1$, we could use this to solve

for any $p_{\mathit{i}}$ to obtain a ranking score for all teams. Instead we will use the point system.

2.6 Template model builder (TMB, R library)

Template model builder is an R/Cpp-library which greatly simplifies the means of estimating hierarchical models. 272829. It makes use of three key concepts: automatic differentiation, laplace method and the (generalized) delta method.

Understanding these concepts are not necessary to make use of the program, but they can be helpful.

2.6.1 Automatic differentiation and dual numbers

Dual numbers are 2-dimensional vectors with the unit vectors $1$ and $\epsilon$ denoted by $(a,b)$ or $a+b\epsilon$, with the added property that the square of the second component is identified with $0$. 30, page 41-43

This has applications for differentiation, where it can be used to find derivatives without differentiating. We define two dual vectors⁸ $\mathbf{u}=(u,u^{\prime })$ and $\mathbf{v}=(v,v^{\prime })$.

We note the similarity between the second component and the rules from differentiation. While the first four rules are trivial, the last requires an explanation. This is found by the tangent expansion⁹ of the function:

The chain rule is also applicable:

By mapping a variable to $x\to (x_0,1)$ and a constant to $c\to (c,0)$, any¹⁰ function can be differentiated by applying the chain rule until sufficiently elementary functions can be computed in order.

For a thorough introduction to the automatic differentiation, we refer to the paper of the stan math library. 31 However, this is only useful for understanding the underlying math of TMB, not for using TMB, so it can safely be ignored.

2.6.2 Laplace method

The second tool is the laplace method, which helps us approximate the marginal distribution and estimate the mean of the fixed parameters and random effects.

figure 5: Normal approximation

The laplace method is based off the tangent series and a quadratic-exponential integral identity¹¹.

We also assume that $f(\theta ,u)$ achieves its maximum at $\widehat{u}$, i.e. ${\partial }_{u}f(\theta ,\widehat{u})=0$, and that $a=b=\infty$ (or that the function decays sufficiently fast from $\widehat{\theta }$). And last, we assume that it achieves it's peak at $\widehat{u}$, i.e. ${\partial }_{uu}^{2}f(\theta ,\widehat{u})<0$.

The proof then follows.

We will use the special case of $M=-1$, thus $f(\theta ,u)=-log(L(\theta ,u))$ and

The multivariate form is slightly different. 32, equation 4TMB uses this to approximate the posterior distribution.

A related result is the posterior central limit theorem¹². We leave out the details and the assumptions, but state the result, often attributed to Bernstein and Von Mises: 33 34

where $\mathcal{I}(\widehat{u})=\bar{n}{\partial }_{uu}^{2}log(p(u|\theta ))$ is the observed information¹³ (i.e. $E(\mathcal{I})=\mathcal{I}$, where $\mathcal{I}$ is called the information) and $n$ is the number of observations.

2.6.3 Delta method

The third tool is the (generalized) delta method,35, page 240-243 which helps us estimate the standard deviation of the fixed parameters and random effects.

The regular delta method states that if there exists a sequence of random variables $X_n$ such that

where $\stackrel{D}{\to }$ denotes convergence in distribution, then for a differentiable function $g$, then

The method used in TMB is a variant that approximates the distribution using the laplace method.

If the posterior distribution of $\lambda |y$ is asymptotically normal with mode $\tilde{\lambda }$, then

where $n$ is the number of observations and $\mathrm{Var}(\lambda )\approx \Sigma =\overline{-{\partial }_{\lambda }^{2}log(p(\lambda =\tilde{\lambda }|y))}$ is the inverse of the negative hessian of the log posterior. Instead of the posterior, one may substitute this for $L(\lambda |y)\times p(\lambda )$, as the factor $p(y)$ is constant (and thus its derivative zero). 36

The more general form where $g(\theta ,u)$ also depends on the random effects $u$, TMB uses a more general estimate:

where $H={\partial }_{uu}^{2}f(\widehat{u}(\theta ),\theta )$, the random effects part of the hessian of the objective function, and $J=D_{\theta }(\widehat{u}(\theta ),\theta )$, the jacobian of $(u(\theta ),\theta )$ wrt. $\theta$. 37

2.6.4 Bugs

TMB does have some quirks and unexpected behaviour. We list some of those here:

• Bad naming: The negative log-distributions are just called distributions. So evaluating $\mathrm{MVNORM}$ yields the negative logarithm of a normal distribution, not the normal distribution.
• Mismatched interfaces: The vectors and matrixes are based off the Eigen library 38, but the arrays are custom made for TMB, and have an incomplete interface compared to vectors and matrixes.
• Unused data and parameter values are silently ignored.

These are practically non-issues, which have simple workarounds, but which one should be aware of. TMB is mature enough to have practical applications, as we demonstrate.

2.7 Attack and defence model

The model we will be using will be a two-stage empirical poisson-log-prior hierarchical distribution. Read that twice. We will be using different distributions on the priors, but they will look similar.

The hierarchical model is shown below

• $\mathit{i}$ and $\mathit{j}$ are team indexes.
• $t$ is the round index.
• $y_{\mathit{i}\mathit{j}t}$ is the number of goals team $\mathit{i}$ scores against $\mathit{j}$ in round $\mathit{j}$.
• ${\lambda }_{\mathit{i}\mathit{j}t}$ is the score parameter for the score.
• $\mu$ is the log-average number of goals overall.
• $\gamma$ is the home advantage. $\mu$ and $\gamma$ are also called fixed effects.¹⁴
• $\alpha$ and $\beta$ are the attack and defence parameters, respectively. We also refer to these as random effects or latent random variables.¹⁵ They attempt to model how each team perform against each other.

Before we describe the priors for the different models, it's important to note the strengths and weaknesses of this model, so we can set realistic expectations of the model performance.

• The model is simple and intuitive to understand.
• It captures the performance of each team individually for every match.
• There isn't a trend parameter; one might expect a team to improve over a season.
• The average and advantage parameters are shared between all teams over all seasons. One would prerably wnat to model each team/season with their separate parameters.

Some of these issues may be resolved by modifying the model, but with a limited data set, and no way to produce more, it's important to keep the model simple to prevent too much overfitting.

Now, for the priors, we make the assumption that these are discretely $\mathrm{Var}(1)$ distributed, or continuously $\mathrm{VAR}$ distributed.

2.7.1 Time independent model

The simplest model is letting the team parameters be constant throughout the season. We will not study this model in detail, but we mention it.

where $w_{\mathit{i}}\sim N(0,\Sigma )$.

2.7.2 Discrete models

In the discrete case we can write the general function

where the absolute value of the eigenvalues of $\Phi$ are smaller or equal than $1$, and $w_t\sim \mathrm{WN}(0,\Sigma )$. In the initial case (unconditional case), we have

where $w_0^{*}=\mathcal{N}(0,{\Sigma }^{*})$, with ${\Sigma }^{*}=\overline{\mathrm{vec}}(\overline{I-{\varphi }_1\otimes {\varphi }_1}\mathrm{vec}({\Sigma }_w))$ in equation (28).

We consider three cases for conditions of $\varphi$:

• $\Phi =0$: The terms in the sequence are independent. This is known as a white noise process, or more specifically a gaussian white noise, as $\epsilon$ are normal distributed. SO this case degenerates to a multivariate normal distribution.
• $\Phi =1$: The next term is the previous term plus a random step. This is known as a random walk.
• Unrestricted. This is a general VAR(1) model. We consider it to be stable as long as the eigenvalues of $\Phi$ are strictly smaller than $1$.

2.7.3 Continuous models

In the continuous case we can write the general function

where the absolute value of the eigenvalues of $\theta$ should be strictly positive in the real part to be stable. In the initial case (unconditional case), we have

where the random integral has the distribution of $\mathcal{N}(0,{\Sigma }^{*})$, with ${\Sigma }^{*}=\overline{\mathrm{vec}}(\overline{\theta \oplus \theta }\mathrm{vec}({\sigma }_w{\sigma }_w^{\top }))$ as in equation (38).

We consider two cases for conditions of $\theta$:

• $\theta =0$: The next term is the previous term plus a random step of a given length. This is known as a wiener process, which is a continuous extension of a random walk, so we refer to this as a continuous random walk.
• Unrestricted. This is a general VAR(1) model. We consider it to be stable as long as the eigenvalues of $\theta$ should be strictly positive in the real part to be stable.
• $\theta \to \infty$ degenerates to a white noise process as in the discrete case, so we ignore this one.

2.8 Result score

2.8.1 Match

While the model itself describse the distribution of the score of a single team, that alone won't help us decide the match winner. For each match we sample from two poisson distributions, or from one skellam distribution.

As usual, the winner is the team with the highest score. If the score of each team is drawn from two poisson distributions, we get $y_{t\mathit{i}}$ and $y_{t\mathit{j}}$. It is then a simple matter of comparing the two, i.e. check which condition holds in $y_{t\mathit{i}}⋛y_{t\mathit{j}}$.

The equivalent condition in terms of the skellam distribution, is to define ${\mathit{k}}_{t\mathit{i}\mathit{j}}=y_{t\mathit{i}}-y_{t\mathit{j}}$ and check ${\mathit{k}}_{t\mathit{i}\mathit{j}}⋛0$. So either of these may be used to determine the winner.

Winning a match gives the winning team three points, and the loser none. A tie gives each team one point. This system is known as three points for a win, and is common in football.

2.8.2 Season

Each team plays against every other team twice, in a double round robin system. If there are $n$ teams, then each team plays $2(n-1)$ matches. The score for each match is added up to a final score, from which the overall seasonal winner is determined.

2.9 Likelihood function

After setting up the model, we want to find the optimal parameters

In the independent case the product of the priors reduce to $\prod _{\mathit{j}=0}^{m}P({\alpha }_{\mathit{j}},{\beta }_{\mathit{j}}|\Sigma )$.

With the accompanying log-likelihood:

By using TMB to apply the laplace approximation to the log-likelihood of the model, we obtain a function of $(\mu ,\sigma ,\varphi ,\Sigma )$ (or $\theta$ in the continuous case), which we can maximize¹⁶. By maximization we obtain $(\widehat{\mu },\widehat{\sigma },\widehat{\varphi },\widehat{\Sigma })$, the arguments maximizing the likelihood, and $(\widehat{\alpha },\widehat{\beta })$, the mode of the posterior, which we use to find the expected value for ${\lambda }_{\mathit{i}}$ for team $\mathit{i}$.

2.10 Measuring model quality

2.10.1 Average score parameter

The average of the score parameters for the matches given time-independent parameters is given by

with $E(\mathrm{avg}({\lambda }_{\mathit{j}}))\approx e^{\mu }$ and $\mathrm{Var}(\mathrm{avg}({\lambda }_{\mathit{j}}))\approx exp(2·0)({\Sigma }_{11}+{\Sigma }_{22})=\mathrm{trace}(\Sigma )$. We remark that $E(f(X))\ne f(E(X))$, so this is only an approximation, nor does it account for dependencies.

This value is not very interesting, as it doesn't help us determine the better team; all lamda-values are equal here. To get comparable lamdas, we instead look at ${\lambda }_{\mathit{j}}|y$. They don't have any nice looking expressions, so we instead use numeric methods to approximate the mean and the variance.

2.10.2 Model selection and fitting

Underfitting means to have a model that isn't sufficiently complex to model the target, i.e. the assumed model is too simple to accurately describe the target.

Overfitting is the opposite: having a model too complex for the target. This often tends to model the noise instead of the underlying distribution.

In both cases the fitted models fail to predict new data points. The aim in model selection is to find an optimal model that neither too strict or too flexible.

In figure 6 we have an example of a target distribution that is modelled by three different polynomials of different orders.

figure 6: Regression fitting

Visually, we can see that the dashed blue model is "best". To find this mathematically we use information criterion values, such as AIC.

2.10.3 Akaike information criterion

AIC is a goodness-of-fit value, giving a lower value for better models. 39 A related value, which we call AIC star, is shown in the below equation:

For small samples, one may subtract a correction term $\frac{\mathit{k}(\mathit{k}+1)}{n-\mathit{k}-1}$ to obtain the corrected criterion ${\mathrm{AICc}}^{*}$, but as this is approximately zero for large samples (assuming $\mathit{k}$ is small), we may ignore this.

The theoretical derivation of AIC includes the quantity $-2log(\widehat{L})$, known as the deviance. So the usual definition of AIC is $\mathrm{AIC}=-2·{\mathrm{AIC}}^{*}$. Because the right hand side of the expression becomes easier to interpret, and it doesn't affect the ranking of the models, we omit the factor $-2$.

For the model in figure 6, we have information criterions for the three fitted models, and for five other polynomial models in figure 7.

figure 7: Information criterion

Using the ${\mathrm{AICc}}^{*}$, we would correctly select the model with the same order as the target, but the coefficients would be different. The $\mathrm{AIC}$ would choose a fourth order curve, which is also close to the true order. Simply relying on the likelihood would have selected an overfitted model; this is the reason for introducing penalizing terms.

The usefulness of $\mathrm{AIC}$ comes from its simple assumptions. It doesn't assume anything about the model; as long as we know the likelihood and the number of parameters, we can calculate the criterion value.

2.10.4 Numerical simulation and estimation of parameters

Using the normal posterior distribution for the attack and defence parameters, we can simulate new values by drawing from the distributions. This can be used to numerically estimate the mean and deviations of the $\lambda$s, and the scores.

From repeated sampling of a season's result, we may obtain a distribution for a team's score for the given attack/defence parameters.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

def season_result(A, gamma, mu):
    scores = np.zeros(shape=teams)

    result_point = {
        -1 : 0,
         0 : 1,
         1 : 3,
    }

    for home, away, hround, around in matches:
        hlamda = np.exp(A[hround,home,0] - A[around,away,1] + gamma + mu)
        alamda = np.exp(A[around,away,0] - A[hround,home,1] - gamma + mu)

        hscore = np.random.poisson(lam=hlamda)
        ascore = np.random.poisson(lam=alamda)

        result = np.sign(hscore - ascore)

        scores[home] += result_point[ result]
        scores[away] += result_point[-result]

    return scores

2.11 Prediction

2.11.1 Naive prediction

While not unexpected, home teams usually have an advantage in matches, often referred to as the home advantage. We model this by the $\gamma$ parameter in our models. Statistically, the home team wins roughly 45 % of the matches in football. 40 In comparison, the home team won 47 % of matches in Eliteserien 2019. 41, so without any knowledge about the teams, a good strategy would be to just bet on the home team. This will be correct almost half the time.

2.11.2 Most likely result outcome

The match result is either a win, tie or a loss.

Using the estimated values ${\lambda }_{\mathit{i}}$ for each team in a given match, we can calculate the loss probability $P_L=P(\mathit{k}>0|{\lambda }_{\mathit{i}1},{\lambda }_{\mathit{i}2})$, the tie probability $P_T=P(\mathit{k}=0|{\lambda }_{\mathit{i}1},{\lambda }_{\mathit{i}2})$, and the win probability $P_W=P(\mathit{k}<0|{\lambda }_{\mathit{i}1},{\lambda }_{\mathit{i}2})$, using the skellam distribution. We then select the most likely result as our guess.

This method has a flaw in that $P_T$ will always be smaller than either $P_L$ or $P_W$, as long as ${\lambda }_{\mathit{i}1},{\lambda }_{\mathit{i}2}\ge 1$ so it will never guess a tie. The proof is simple: $P_T({\lambda }_{\mathit{i}1}=1,{\lambda }_{\mathit{i}2}=1)<\frac{1}{3}$, and is a maximum. I.e. increasing either ${\lambda }_{\mathit{i}1}$ or ${\lambda }_{\mathit{i}2}$ will make this value smaller. So either $P_L$ or $P_W$ must be greater than $\frac{1}{3}$, and be our guess.

However, around a third of matches result in a tie, 42, but figure 2 shows the average of ${\lambda }_{\mathit{i}}$ to be around $1.5$, so unless the model is really accurate and can estimates below $1$, this method will be wrong approximately one third of the time.

2.11.3 Most likely score outcome

The match score is a pair of number, e.g. $1-4$, $2-2$ or $7-1$.

To make tie guesses more likely, we may want to use the mode instead. The mode represents the most likely score outcome, so instead of looking at what's most likely of a win, tie or loss, we look at each score outcome individually.

This makes sense because we fit the model to score, not the result of a match. However, it should more heavily favour ties, even when winning results combined would be more likely.

2.11.4 Weighted outcome

Another method that tries to correct the win bias, is to apply weights to the win, tie and loss probabilities. We then select the result with the highest weighted probability.

3 Data analysis

3.1 Eliteserien 2019

The resulting scores from 2019 are taken from 43. The scores are displayed in table 1.