GAUSS-NEWTON algorithm to solve a non-linear system?

Question

I need a parameter to use in Poisson. I want to calculate this parameter using the Gauss-Newton algorithm, starting from a nonlinear system. I want to point out that i am assuming that the distribution is Poisson with parameter $\mu$

My goal is to use Poisson (which is not relevant to this question, because with this question i just want to get the only parameter of Gauss-Netwon to use in Poisson) to individually calculate the probabilities of a football club scoring 0 goals, 1 goal, 2 goals, 3 goals, 4 goals, etc. etc. So I need to find a parameter using the Gauss-Newton algorithm.

In an experimental test, i tried to use the Arithmetic Mean in Poisson, then i use the parameter 1.6 to calculate the probability that a football team individually scores goals, 1 goal, 2 goals, 3 goals, 4 goals, etc. etc. The results are decent, but I would like to get something more precise. So I thought of using another parameter, which is what i would like to find with the Gauss-Newton algorithm and the nonlinear system.

The basic, initial idea, is to start from List_Goals: [1, 2, 2, 1, 2] which i organize in the non-linear system, building it with the goals [1, 2, 2, 1, 2] and the number of times the event occurs

BASIC DATA:

List_Goals: [1, 2, 2, 1, 2]
Matches played: 5 
Standard Deviation: 0.55

NON-LINEAR SYSTEM:

If a football club, in the last 5 matches,scores [1, 2, 2, 1, 2], it means that:

• 0 goals, scored 0 times ( );
• 1 or less goals, scored 2 times (1, 1);
• 2 or less goals, scored 5 times (1, 2, 2, 1, 2);

therefore:

If the goals scored are <= 0 events, then i will have 0/5 = 0;
If the goals scored are <= 1 events, then i will have 2/5 = 0.4;
If the goals scored are <= 2 events, then i will have 5/5 = 1;

f(0) = 0/5 = 0;
f(1) = 2/5= 0.4;
f(2) = 5/5= 1;

System: {f(0) = F(0) } therefore 0/5= 0
        {f(1) = F(1) } therefore 2/5= 0.4
        {f(2) = F(2) } therefore 5/5= 1

In this case i have three equations in one unknown, but obviously the system must be set up with a number of equations ranging from F(0) to F(max_goals_scored). This way I would have max_goals_scored + 1 equations.I should therefore start the solution starting from max_goals_scored+1 equations and increase, for example up to a maximum of 5 goals:

F(0), F(1), F(2), F(3), F(4), F(5) 

GAUSS-NEWTON ALGORITHM:

I'm having trouble using the Gauss-Newton algorithm, so i'm asking for help. I understand the theory, but i'm having trouble putting it into practice. Specifying that, what I'm looking for for the resolution and the final result is this:

a) I identify the Jacobian matrix J, which in general is composed of the partial derivatives with respect to all the parameters. In this case i only have one parameter and the matrix is actually a column vector;

b) I calculate the transpose of J multiplied by J, which in this case is a scalar (formally a number, which however in this case is a function of the unknown parameter);

c) algorithmically speaking it looks good: the inverse of J^T*J is still a scalar, so I have to divide by that number (which however will contain the mu parameter);

d) I write the iterative method according to the formula and a stopping criterion;

Naturally, these steps therefore, include functions and functional operators (the derivative). How can i use the Gauss-Network algorithm on the non-linear system and get the final result? Can you show me the various steps with my numerical data applied? (and not just with the formulas). Thank you all!!!

Answer 1

This is not an answer, but is too long for a comment.

It is not entirely clear to me what you are asking, you seem to have all the information necessary. However the notation could be improved.

I am using $F_k(\mu) = e^{-\mu} \sum_{j=0}^k { \mu^j \over j!}$ and am assuming that $f_e(k)$ represents the corresponding experimental datum. It is straightforward to compute $F_k'(\mu)$.

I am guessing that you have a system of equations $r:\mathbb{R} \to \mathbb{R}^3$, $r_k(\mu) = F_k(\mu)-f_e(k)$, then you would like to solve $r(\mu) \approx 0$.

To apply at iteration of Gauss Newton at some iteration $e_n$ you linearise $r$ as in $r(\mu_{n+1}) \approx r(\mu_n) + Dr(\mu_n) (\mu_{n+1}-\mu_n)$ and (attempt to) solve the system $r(\mu_n) + Df(\mu_n) (\mu_{n+1}-\mu_n) = 0$ to compute $e_{n+1}$.

Using ordinary least squares, we have the update rule $\mu_{n+1} = \mu_n - (Dr(\mu_n)^T Dr(\mu_n))^{-1} Dr(\mu_n)^T r(\mu_n)$.

In this case (assuming that I understand correctly) you have, $Dr(\mu) = \begin{bmatrix} F_0'(\mu) \\ F_1'(\mu) \\ F_2'(\mu) \end{bmatrix}$.

A readable explanation of Gauss Newton can be found here, some standard improvements (Armijo step size, LevenbergMarquardt) can be found here.

Another approach:

The following plots the sum of errors squared vs. the Poisson $\mu$ parameter for the data [1, 2, 2, 1, 2]. There could easily be errors in my implementation.

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
import math

goals = [1, 2, 2, 1, 2]
def f(k):
    """Observed data."""
    return len(list(filter(lambda j: j <= k, goals)))/len(goals)

def poisson_cdf(mu, k):
    """Poisson CDF."""
    return sum(mu**j/math.factorial(j) for j in range(0, k+1))*math.exp(-mu)

def lsq_error(mu, k_vals):
    """Sum of squares of errors."""
    return sum((poisson_cdf(mu, k)-f(k))**2 for k in k_vals)

# Find the error for values of k in k_vals.
k_vals = [0, 1, 2]
# Plot x range.
mu_vals = np.arange(0.0, 5.0, 0.1)
lsq_error_vals = np.array(list(lsq_error(mu, k_vals) for mu in mu_vals))

fig, ax = plt.subplots()
ax.plot(mu_vals, lsq_error_vals)

ax.set(xlabel='mu', ylabel='lsq error',
       title='sum of errors squared vs. mu')
ax.grid()

fig.savefig("error.png")
plt.show()

This gives the following plot from which we can visually see that there is a $\min$ around $1.6$.