# Problem with data optimization after sum of squared errors in Python

I have a script that allows you to optimize some data (goals from all matches of a football team), after plotting the sum of the squared errors. The aim is to find a final parameter that is more precise than the arithmetic mean, because i will later use this parameter in the Poisson distribution to calculate the specific and individual numbers of exact goals. Initially I used the Poisson average, but I would like a more precise and better parameter.

I'm not sure if the script is correct. The problem is that I'm not sure if the final parameter is correct (1.69), because i'm not sure that the step of plotting the sum of squared errors with respect to the Poisson parameter μ is what i really need or if I need something different (maybe I should have solved the least squares problem?), OR i'm also not sure if the optimization with BFGS is what I really need.

For example, considering the list goals = [1, 2, 2, 2, 1] which has an average of 1.6, then the final parameter of the script is 1.69.

What do you think is wrong with my code and how could I improve it?

I divide the code into two parts. The first part plots the sum of squared errors versus the μ Poisson parameter for the given goals. The second part uses a BFGS local optimizer. I preferred not to enter the limits (bounds) because the minimum and maximum values vary continuously and are therefore always unknown (goals = [1, 2, 2, 2, 1] is only a example). I would have liked to use an optimizer/solver like Newton-Gauss, but I was advised to use any solver because I am solving a problem with only one dimension and not hundreds.

import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import poisson
from scipy.optimize import minimize

goals = [1, 2, 2, 2, 1]

# 1 PART
# Plots the sum of squared errors versus the μ Poisson parameter for the given goals
def f(k):
"""Observed data."""
observed_data = len(list(filter(lambda j: j <= k, goals)))/len(goals)
return observed_data

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.
max_value = max(goals)
k_vals= list(range(max_value + 1))
print("Number of Functions: ", k_vals)

###############################################################
# 2 PART
# optimize with BFGS local optimizer
result = minimize(lambda x: lsq_error(x[0], k_vals), x0=[0])
print("Final parameter: ", result.x[0]) #1.6931671442067795


Thank you all!