I was seeking a formula to relate goals scored and goals allowed to points in association football.

In association football, three points are awarded to the team winning a match, with no points awarded to the losing team. If the game is drawn, each team receives one point.

Firstly, I worked with points percentage — points divided by available points, which I set to the number of matches multiplied by adjusted 3-pointer. The general form was: $${\rm points}\,{\rm percentage}= \frac{{\rm goals}\,{\rm scored}^{\mathcal{A}}}{{\rm goals}\,{\rm scored}^{\mathcal{A}}+ {\rm goals}\,{\rm allowed}^{\mathcal{B}}}$$ The league tables' total points percentage is almost sure not an integer.
This is Pythagorean expectation (only for wins/losses, not for draws) with two forced different exponents! Football presents a Pythagorean challenge because about a quarter of its matches end in draws (*related to trimean ? no.), worth one point each in the league tables. Because wins are worth three points, the most accurate Pythagorean analyses overestimate most teams’ points totals, as football analyst Howard Hamilton has found.
I call adjusted 3-pointer instead of three because how many average points each team receives are not same. I want to achieve 'reasonable' relationship among adjusted 3-pointer (not parameter, but formula like Pythagenport, Pythagenpat), goals scored and goals allowed for each team: $${\rm points}_{predicted}= {\rm points}\,{\rm percentage}\cdot{\rm adjusted}\,3{\rm -pointer}$$ If you treat adjusted 3-pointer as parameter:
Here is the win estimator (Python Code)

# 17 Th12 '12

import scipy.optimize
import numpy as np
z=np.array([36,42,29,24,27,29,23,27,24,23,22,20,25,16,17,15,18, 9,15,10])
pred=lambda w: (x**w[0]/(x**w[0]+y**w[1]))*(2+1/(w[2]**2+1))*N
w,_=scipy.optimize.leastsq(lambda w: (z-pred(w)),(1,1,.0009332150))

Hence, we obtain

adjusted 3-pointer=2.742957368624336

As you can see, I choose ${\rm adjusted}\,3{\rm -pointer}= 2+ 1/(w^{2}+ 1)\in\left ( 2, 3 \right )$ because the average total points of football matches is on this interval. And ${\rm adjusted}\,3{\rm -pointer}\lessapprox 2.75$ is scientific* (noted above). On the other side, $w> {\rm Best}\,w\Rightarrow 2\lessapprox{\rm adjusted}\,3{\rm -pointer}$ leads to many under-performing teams, despite of well RMSE, I still want adjusted 3-pointer as a formula (function), but parameter. I have tried to turning that into an ODE (In fact, a PDE.) but unsuccessfully. I need your help to find it.


1 Answer 1


In order to apply ordinary least squares method in your case and to try and answer your questions afterwards, I first derive two types of equations:

1. Linear observation equations. Denote the expected points percentage by $p$, total goals scored by $TGS$, total goals allowed by $TGA$, and the corresponding exponents by $\alpha$, $\beta$ respectively. Then your first formula becomes \begin{equation} p = \frac{\text{TGS}^{\alpha}}{\text{TGS}^{\alpha} + \text{TGA}^{\beta}} \tag{1} \end{equation}

We can derive a model that is linear in exponents $\alpha, \beta$ in the following way. First, notice that $$ \frac{p}{1-p} = \frac{\text{TGS}^{\alpha}}{\text{TGA}^{\beta}} $$ Taking the logarithm $$ \log{\left(\frac{p}{1-p}\right)} = \alpha \cdot \log{\left(\text{TGS}\right)} - \beta \cdot \log{\left(\text{TGA}\right)} $$ Now, if we introduce $$ y := \log{\left(\frac{p}{1-p}\right)} \qquad x_1 := \log{\left(\text{TGS}\right)} \qquad x_2 := \log{\left(\text{TGA}\right)} $$ we arrive at a simple linear relationship: $$ y = \alpha \cdot x_1 - \beta \cdot x_2 $$ The dependent variable $y$ is called a response and independent variables $x_1, x_2$ are known as regressors.

Given some data $y_{i}$, $x_{i1}$, $x_{i2}$, $i = 1, \dots n$, we can write a linear model in parameters $\alpha, \beta$ as \begin{equation} y_i = \alpha \cdot x_{i1} - \beta \cdot x_{i2} + \varepsilon_{i} \tag{2} \end{equation}

Under assumptions that $\varepsilon$ has mean $0$ and $\varepsilon \perp x_{1}, x_{2}$, we can use ordinary least squares (OLS) in order to find the estimates $\widehat{\alpha}$ and $\widehat{\beta}$. This gives us the estimate of $y_{i}$ as $$ \widehat{y_{i}} = \widehat{\alpha} \cdot x_{i1} - \widehat{\beta} \cdot x_{i2} $$

We can then also express the estimate for $p_{i}$ as \begin{equation} \widehat{p_{i}} = \frac{1}{1 + \exp{\left(-\widehat{y_{i}}\right)}} \end{equation}

Finally, the mean squared error (MSE) of an estimator of $p$ is defined as \begin{equation} \text{MSE} = \frac{1}{n} \sum_{i = 1}^{n} \left( p_{i} - \widehat{p_{i}} \right)^{2} \tag{3} \end{equation} And taking the square root of the MSE yields the root mean squared error (RMSE).

2. Nonlinear observation equations. If we introduce a new parameter $\gamma$ and define adjusted 3-pointer as a real valued function $h(\gamma)$. For example, as you defined in your question \begin{equation} h(\gamma) = 2 + \frac{1}{1 + \gamma^{2}} \tag{4} \end{equation} then the formula (1) becomes: \begin{equation} p = \frac{TGS^{\alpha}}{TGS^{\alpha} + TGA^{\beta}} \cdot h(\gamma) \tag{5} \end{equation}

The response $y = \log{\left( p / (1 - p) \right)}$ is now some non-linear function $f: \mathbb{R}^{3} \rightarrow \mathbb{R}$ of parameters $\pmb{w} := [\alpha, \beta, \gamma]^{T}$. We can again derive a model that is linear in parameters $\pmb{w}$ in the following way.

Let $\Delta y = f(\pmb{w}) - f(\pmb{w_0})$ and $\nabla f (\pmb{w_{0}})$ a vector of partial derivatives $\left[ \frac{\partial f}{\partial \alpha}, \frac{\partial f}{\partial \beta}, \frac{\partial f}{\partial \gamma} \right]^T$ evaluated at $\pmb{w_0} = [\alpha_0, \beta_0, \gamma_0]^T$. Linearization of $y = f(\pmb{w})$ yields the linear model in $\Delta \pmb{w} = \pmb{w} - \pmb{w_0}$: \begin{equation} \Delta y = \nabla f (\pmb{w_{0}})^T \Delta \pmb{w} + \varepsilon \tag{6} \end{equation}

Again, under assumptions that $\varepsilon$ has mean $0$ and $\pmb{\varepsilon} \perp \pmb{w}$, we can use ordinary least squares (OLS) to find the estimate $\widehat{ \Delta \pmb{w}}$, and use it to get $\widehat{\pmb{w}} = \pmb{w_0} + \widehat{ \Delta \pmb{w}}$.

Iterative scheme. To get a final estimate $\widehat{\pmb{w}}$, start at some initial guess $\pmb{w_0} = [\alpha_0, \beta_0, \gamma_0]^T$. Solve (6) using OLS to get an update $\Delta\widehat{\pmb{w}}$ and use it to get the new estimate: $$ \widehat{\pmb{w_{1}}} = \pmb{w_{0}} + \Delta\widehat{\pmb{w}} $$ Set $\pmb{w_{0}} \leftarrow \pmb{w_{1}}$ and repeat, until $\newcommand{\norm}[1] {\left\lVert#1\right\rVert}\norm{\widehat{ \Delta \pmb{w} }}_{2}$ becomes sufficiently small.

Now, to try and answer your questions. Yes, using a model based on nonlinear observation equations derived from (5), can lead to an estimator $\widehat{p}$ with a better RMSE, since a more complex model with more parameters can better fit the data. But, I don't understand what you mean by:

On the other side, $w> {\rm Best}\,w\Rightarrow 2\lessapprox{\rm adjusted}\,3{\rm -pointer}$ leads to many under-performing teams

This can mean that $\widehat{p}$ is a biased estimator of the expected points percentage $p$. To get an unbiased estimator $\widehat{p}$ using OLS, the assumptions of OLS must be satisfied. Check that the estimated residuals $\widehat{\varepsilon}_i = y_i - \widehat{y_i}$ have mean $0$ and are not correlated with regressors $x_i$. Otherwise the iterative scheme of updates $\Delta\widehat{\pmb{w}}$ is not valid. When using some package, e.g. scipy.optimize, check if leastsq returns any warnings. Or you can try implementing the iterative scheme yourself and check the assumptions. If assumptions are not satisfied, you have to choose a different model. Also, a more complex model can easily overfit the data, so cross-validation is usually used to validate the results.

I have tried to turning that into an ODE (In fact, a PDE.) but unsuccessfully. I need your help to find it.

To model the points percentage $p$ with a PDE, you need to describe the relation between partial derivatives of $p$ with respect to the parameters $\alpha, \beta, \gamma$. From your description of $p$, I don't see that it has to satisfy some PDE.

Final remarks.

  • If you are looking for this relationship for $p$ in a data driven way, what your question points to, this is a standard problem in machine learning. But the result will not be some simple formula, but a fitted model.
  • For example, you can fit a neural network to find the relationship you are looking for. Since even a single hidden layer network can approximate any continuous function. Again, to not overfit, validate the result, for example by using cross-validation.
  • But, when given little amount of data, avoid fitting a more complex model. To model rare events such as goals, the Poisson distribution is typically used.
  • Finding an analytical expression that matches data from an unknown function is in principle NP-hard. Instead, any polynomial and rational function modeling is typically knowledge-driven.
  • This means that you have to come up with a simple formula for the relationship you are looking for based on the nature of the problem and then find the values of parameters using some optimization technique.
  • $\begingroup$ What style! This is the answer I have waited for a long long time. Thanks a real lot! $\endgroup$
    – user822157
    Commented Sep 7, 2022 at 1:07
  • 1
    $\begingroup$ Thank you, I am happy to hear that! $\endgroup$ Commented Sep 7, 2022 at 20:05
  • $\begingroup$ Would you like to help me_ math.stackexchange.com/q/4502888/822157 ? I have protected it for a long time. $\endgroup$
    – user822157
    Commented Sep 25, 2022 at 10:28

You must log in to answer this question.