# Solve a non-linear system, in Python, with the GAUSS-NEWTON algorithm? (Jacobian matrix J, etc.)

I would like to solve a non-linear system (which contains the goals of a football team in previous matches) using the Gauss-Netwon algorithm, in order to find the parameter (of frequency) to use as lambda in the simple Poisson Distribution to calculate the probability of scoring a given number of goals. The question is not about Poisson, but only about the Gauss-Newton algorithm. I would need the executable solution as Python code, because with just a few suggestions i wouldn't be able to solve the problem. Is there any library that uses the Gauss-Newton algorithm or do you have to write everything manually? Can you show me how to solve the nonlinear system with Gauss-Newton with Poisson please?

I have a nonlinear system with three equations in a single unknown. The nonlinear system was created using the analytical expression of Poisson CDF with the Distribution Function. I calculate the cumulative frequencies from the observations and use them as data. The system works properly, there are no problems. In the system are the goals of a football team in previous matches.

F(t) := P(X <= t) ~ sum_i_frequency(observation_i <=t) / total_observation =: f(t)

List_Goals: [1, 2, 2, 1, 2]
Matches played: 5

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;


The result of this system is: 0, 0.4, 1

The system runs successfully with Python code:

import numpy as np
Goals = [1, 2, 2, 1, 2]
probability_mass_function = np.bincount(Goals)/len(Goals)

cummulative_mass_function = probability_mass_function.cumsum()

print("probability_mass_function: ", probability_mass_function)
#result: ([0. , 0.4, 0.6])

print("cummulative_mass_function: ", cummulative_mass_function)
#result: ([0. , 0.4, 1. ])


## GAUSS-NEWTON

Now, i would like to solve the nonlinear system using the Gauss-Netwon algorithm, in order to find the parameter to use as lambda in the simple. Poisson distribution. So I would like to solve the nonlinear system using the following Gauss-Newton algorithm:

So what I would need is:

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;

So I get the parameter (of frequency) to use in the Poisson distribution to calculate 0 goals, 1 goal, 2 goals, 3 goals, etc. ​ Naturally, these steps therefore, include functions and functional operators (the derivative). How can I use the Gauss-Network algorithm (with Python) on the non-linear system and get the final parameter I can use in Poisson distribution?

I would need the executable solution as Python code, because with just a few suggestions I wouldn't be able to solve the problem.

Thank you all!

Python mathematical libraries numpy and scipy have routines for solving systems of linear equations: numpy.linalg.solve and scipy.linalg.solve.

These are based on Gaussian elimination, rather than invertsing matrix (which can be achieved, e.g., by numpy.linalg.inv.)

See also Solving equations and inverting matrices for the full list of available routines.

• Sorry, I just wanted to warn you that I don't find your answer useful, because it doesn't answer what I'm looking for. As you have seen there is a Bounty +50. If you help me by editing your answer and writing a better one, I will be happy to give you the bounty. Thank you Dec 9, 2023 at 0:28
• FOR READERS: Of course I would also like to read solutions from other users :) Dec 9, 2023 at 0:28

I hope this is what you are looking for

I am using scipy.optimize.least_squares method which implements Gauss-Newton method. (Doc - https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)

An example code I tried is below(Edited after clarification from comment)

import numpy as np
from scipy.optimize import least_squares
from scipy.stats import poisson

# Define the Poisson distribution model
def poisson_model(params, x_data):
# params: frequency parameter for the Poisson distribution
return np.cumsum(poisson.pmf(x_data, params))

# Define the objective function
def objective_function(params, x_data, y_data):
# x_data: the number of goals in previous matches
# y_data: observed cumulative frequencies
# Poisson distribution cumulative probabilities
expected_cumulative_probabilities = poisson_model(params[0], x_data)
# Compute the residuals
residuals = expected_cumulative_probabilities - y_data
return residuals

# Example data: number of goals in each match
goals_data = np.array([1, 2, 2, 1, 2])
max_goals = np.max(goals_data)
goal_counts = np.arange(0, max_goals + 1)
observed_cumulative_frequencies = np.cumsum(np.histogram(goals_data, bins=np.arange(-0.5, max_goals + 1.5))[0])

# Initial guess for the frequency parameter
initial_frequency = 1.0
initial_params = np.array([initial_frequency])

# Use least_squares to perform Gauss-Newton optimization
result = least_squares(objective_function, initial_params, args=(goal_counts, observed_cumulative_frequencies))

# The optimized frequency parameter is in result.x
optimized_frequency = result.x[0]

print("Optimized Frequency Parameter:", optimized_frequency)


In this code, observed_cumulative_frequencies represents the cumulative frequencies of goals. For example, if the observed cumulative frequencies for each goal count are [0, 1, 4, 5, 5, 5], it means that there are no occurrences of less than 0 goals, 1 occurrence of less than 1 goal, 3 occurrences of less than 2 goals, and so on.

• What is observed_frequencies = np.array([0, 1, 3, 1, 0])? Why did you write 0, 1, 3, 1, 0? Can you explain to me with an example? Dec 10, 2023 at 5:55
• Looks like I didn't correct that, sorry,In your post, you have written goals as [1,2,2,1,2]. However you just need unique values so goals will be now [1,2] and then observed frequency will be [2,3] since 1 has a frequency of 2 and 2 has a frequency of 3 Dec 10, 2023 at 6:03
• I'm not understanding. Where should I write my [1, 2, 2, 1, 2]? Dec 10, 2023 at 6:10
• Okay, you have your data as football goals. matches played : 5 and goals in those 5 matches are [1,2,2,1,2]. From this, you have to calculate frequency. In 5 matches, there were 2 matches in which goals scored were 1 and there were 3 matches in which goals scored were 2, so you goal data is now [1,2] and observed frequency is [2,3]. I hope this is clear. Let me know if you need more clarification Dec 10, 2023 at 6:16
• There are two problems.You are counting differently than I would like. In my question I count differently.You counted 2 games in which 1 goal was scored, and then 3 games in which 2 goals were scored. As in my question, I would like to count 0 matches in which 0 or fewer goals were scored; 2 matches in which 1 or fewer goals were scored; 5 matches in which 2 or fewer goals were scored. Also, in your code, I don't have the possibility to DIRECTLY insert my List_Goals: [1, 2, 2, 1, 2].Could you fix your question a little and improve it please? I will be happy to deliver the bounty to you Dec 10, 2023 at 6:26