Решить нелинейную систему на Python с помощью алгоритма ГАУССА-НЬЮТОНА? (Матрица Якоби J и т.д.)

Я хотел бы решить нелинейную систему (которая содержит голы футбольной команды в предыдущих матчах) с использованием алгоритма Гаусса-Ньютона, чтобы найти параметр (частоты), который будет использоваться в качестве лямбда в упрощенном распределении Пуассона для расчета вероятности забить заданное число голов. Вопрос не о распределении Пуассона, а только о методе Гаусса-Ньютона. Мне нужна исполнимая реализация в виде кода на Python, потому что несколько подсказок не помогут мне решить эту задачу. Существует ли библиотека, которая использует алгоритм Гаусса-Ньютона, или его нужно реализовать вручную? Можете показать, как решить нелинейную систему с помощью метода Гаусса-Ньютона и распределения Пуассона, пожалуйста?

У меня есть нелинейная система с тремя уравнениями и одной неизвестной. Нелинейная система была создана с использованием аналитического выражения для функции распределения Пуассона с функцией распределения. Я рассчитываю накопленные частоты на основе наблюдений и использую их как данные. Система работает правильно, проблем нет. В системе учтены голы футбольной команды в предыдущих матчах.

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

Список голов: [1, 2, 2, 1, 2]
Матчей сыграно: 5

Если голов забито <= 0, тогда: 0/5 = 0;
Если голов забито <= 1, тогда: 2/5 = 0.4;
Если голов забито <= 2, тогда: 5/5 = 1;

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

Система: {f(0) = F(0)} таким образом 0/5 = 0;
        {f(1) = F(1)} таким образом 2/5 = 0.4;
        {f(2) = F(2)} таким образом 5/5 = 1;

Результат этой системы: ;0, 0.4, 1

Система успешно выполняется с помощью Python кода:

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)
#результат: ([0. , 0.4, 0.6])

print("cummulative_mass_function: ", cummulative_mass_function)
#результат: ([0. , 0.4, 1. ])

ГАУСС-НЬЮТОН

Теперь я хотел бы решить нелинейную систему с использованием алгоритма Гаусса-Ньютона, чтобы найти параметр для использования в упрощенном. распределении Пуассона. Итак, я хотел бы решить нелинейную систему, используя следующий алгоритм Гаусса-Ньютона:

Итак, что мне нужно:

a) Я определяю Якобиан матрицы J, который в общем случае состоит из частных производных по всем параметрам. В данном случае у меня только один параметр, и матрица на самом деле является вектором-столбцом;

b) Я вычисляю транспонированную J, умноженную на J, что в данном случае является скаляром (формально числом, которое, однако, в этом случае является функцией от неизвестного параметра);

c) алгоритмически говоря, это кажется правильным: обратная от J^T*J все еще является скаляром, поэтому мне нужно разделить на это число (которое, однако, будет содержать параметр mu);

d) Я записываю итерационный метод в соответствии с формулой и критерием остановки;

Таким образом, я получаю параметр (частоты), который можно использовать в распределении Пуассона для расчета 0 голов, 1 гола, 2 голов, 3 голов и т. д.

Естественно, эти шаги, следовательно, включают функции и функциональные операторы (производную). Как я могу использовать алгоритм Гаусса-Ньютона (с Python) на нелинейной системе и получить итоговый параметр, который я могу использовать в распределении Пуассона?

Мне нужна исполнимая реализация в виде кода Python, потому что несколько подсказок не помогут мне решить эту задачу.

Спасибо всем!

Математические библиотеки Python numpy и scipy имеют процедуры для решения систем линейных уравнений: numpy.linalg.solve и scipy.linalg.solve.

Эти методы основаны на методе исключения Гаусса, а не на обращении матриц (что можно сделать, например, с помощью numpy.linalg.inv).

Смотрите также Решение уравнений и обращение матриц для полного списка доступных процедур.

Надеюсь, это то, что вы ищете

Я использую метод scipy.optimize.least_squares, который реализует метод Гаусса-Ньютона. (Документация – https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)

Пример кода, который я попробовал, ниже (исправлен после уточнения из комментария)

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

# Определите модель распределения Пуассона
def poisson_model(params, x_data):
    # params: параметр частоты для распределения Пуассона
    return np.cumsum(poisson.pmf(x_data, params))

# Определите целевую функцию
def objective_function(params, x_data, y_data):
    # x_data: количество голов в предыдущих матчах
    # y_data: наблюдаемые накопленные частоты
    # Кумулятивные вероятности распределения Пуассона
    expected_cumulative_probabilities = poisson_model(params[0], x_data)
    # Вычислите остатки
    residuals = expected_cumulative_probabilities - y_data
    return residuals

# Пример данных: количество голов в каждом матче
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_frequency = 1.0
initial_params = np.array([initial_frequency])

# Используйте least_squares для выполнения оптимизации методом Гаусса-Ньютона
result = least_squares(objective_function, initial_params, args=(goal_counts, observed_cumulative_frequencies))

# Оптимизированный параметр частоты находится в result.x
optimized_frequency = result.x[0]

print("Оптимизированный параметр частоты:", optimized_frequency)

В этом коде observed_cumulative_frequencies представляет собой накопленные частоты голов. Например, если наблюдаемые накопленные частоты для каждого количества голов равны [0, 1, 4, 5, 5, 5], это означает, что нет событий менее 0 голов, 1 событие менее 1 гола, 3 события менее 2 голов и так далее.