Curve fit sum of poisson function to histogram data

Question

I am new here so I hope this post is legible/set up correctly. I am having difficulty trying to use curve_fit in python to fit histogram data to a sum of Poisson distributions. The data is a set of photon counts for different numbers of atoms. The set in this code is truncated to 0 atoms (background) 1 atom or 2 atoms. Whether I use a user defined Poisson or the poisson.pmf function from scipy I am unable to get a successful fit without errors.

```python
from scipy.special import factorial
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np

### User Defined poisson
def poissonian( x, A0, x0, A1, x1, A2, x2 ):
    return (A0 * x0**x * np.exp(-x0) / scipy.special.factorial(x)\
        +A1 * x1**x * np.exp(-x1) / scipy.special.factorial(x)\
        +A2 * x2**x * np.exp(-x2) / scipy.special.factorial(x))

### Built in poisson function
#def poissonian( x, A0, x0, A1, x1, A2, x2 ):
#    return (A0*poisson.pmf(x,x0)+A1*poisson.pmf(x,x1)+A2*poisson.pmf(x,x2))

binwidth = 75
bin2 = 18
xmin = 0
xmax = bin2 * binwidth

counts, bins, bars = plt.hist( counts_77_2x2, bins=bin2, range=(xmin, 
xmax) )
#counts_77_2x2 is a dataset for the number of photons detected from atoms 
#trapped in an optical tweezer. The point of this fit is to determine the 
#probability of loading n atoms in a given optical tweezer, 77 is for tweezer 
#7, 2x2 is the box pixel size for the target tweezer site imaged by emccd


bins = np.rint( bins[1:] - 75 / 2 ).astype( np.int32 )

#for reference
#bins =[  38  112  188  262  338  412  488  562  638  712  788  862  938 
#1012 1088 1162 1238 1312]
#counts= [ 2.  0.  1.  2.  2.  5.  3. 10. 11.  8.  6.  4. 10.  8.  8.  9.  
#7.  3.]


## Initial Guess
p0 = [1,400, 9,650, 8,1050]

#Bounds for fit
bounds=([1,1, 5,525, 5,850], [6,524, 14,800, 14,1350])

coeff, var_matrix = curve_fit(poissonian, bins, counts, p0=p0, 
maxfev=1000000)

print( coeff )```

Using a user defined poisson function yields

....RuntimeWarning: overflow encountered in powerreturn (A0 * x0**x * np.exp(-x0) / scipy.special.factorial(x)....

I have read the first few pages of responses to similar issues, I replaced the factorial with the scipy.special.factorial as I read it would accommodate higher value inputs but it was unable to help.

Using the built in poisson.pmf probability mass function yields

...\anaconda3\lib\site-packages\scipy\optimize\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated warnings.warn('Covariance of the parameters could not be estimated',

To reduce the calculation overhead I scaled down the x values to increments of 1 (all bins are equally spaced so it can scaled up for real peak centers later) Using this I was finally able to get a working (albeit ugly) fit.

```python
import scipy
from scipy.stats import poisson
from scipy.special import factorial
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np


def poissonian( x, A0, x0, A1, x1, A2, x2 ):
    return (A0 * x0**x * np.exp(-x0) / scipy.special.factorial(x)\
        +A1 * x1**x * np.exp(-x1) / scipy.special.factorial(x)\
        +A2 * x2**x * np.exp(-x2) / scipy.special.factorial(x))

binwidth = 75
bin2 = 18
xmin = 0
xmax = bin2 * binwidth

counts, bins, bars = plt.hist( counts_77_2x2, bins=bin2, range=(xmin, xmax) )

#for reference
#bins =[  38  112  188  262  338  412  488  562  638  712  788  862  938 1012 
#1088 1162 1238 1312]
#counts= [ 2.  0.  1.  2.  2.  5.  3. 10. 11.  8.  6.  4. 10.  8.  8.  9.  7.  
#3.]


### set bin separation to one to reduce calculation overhead (bins1)

bins1 = np.linspace( 1, 18, 18 )

## Scaled down guess values
p0 = [1,4, 8,10, 8,15]

#scaled down bounds
bounds=([0.1,1, 0.1,7, 0.1,13], [5,6, 12, 3, 12,18])

coeff, var_matrix = curve_fit(poissonian, bins1, counts, p0=p0, 
maxfev=1000000)

print( coeff )
```

Now that I got a scaled down version working, I wanted to know, does the code using the real data set have an error in it or is there a fundamental limitation for python's ability to fit poissonians for non small x valued data.

Thank you for your time!