Fitting (a Gaussian) With Scipy Vs. Root Et Al

I have now multiple times stumbled upon that fitting in python with scipy.curve_fit is somehow a lot harder than with other tools like e.g. ROOT (https://root.cern.ch/) For example

Solution 1:

Without actual data it is not that easy to reproduce your results but with artificially created noisy data it looks fine to me:

That is the code I am using:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# your gauss functiondefgauss(x, a, x0, sigma):
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

# create some noisy data
xdata = np.linspace(0, 4, 50)
y = gauss(xdata, 2.5, 1.3, 0.5)
y_noise = 0.4 * np.random.normal(size=xdata.size)
ydata = y + y_noise
# plot the noisy data
plt.plot(xdata, ydata, 'bo', label='data')

# do the curve fit using your idea for the initial guess
popt, pcov = curve_fit(gauss, xdata, ydata, p0=[ydata.max(), ydata.mean(), ydata.std()])

# plot the fit as well
plt.plot(xdata, gauss(xdata, *popt), 'r-', label='fit')

plt.show()

As you, I also use p0=[ydata.max(), ydata.mean(), ydata.std()] as an initial guess and that seems to work fine and robust for different noise strengths.

EDIT

I just realized that you actually provide data; then the result looks as follows:

Code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


def gauss(x, a, x0, sigma):
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

ydata = np.array([2., 2., 11., 0., 5., 7., 18., 12., 19., 20., 36., 11., 21., 8., 13., 14., 8., 3., 21., 0., 24., 0., 12.,
0., 8., 11., 18., 0., 9., 21., 17., 21., 28., 36., 51., 36., 47., 69., 78., 73., 52., 81., 96., 71., 92., 70., 84.,72.,
88., 82., 106., 101., 88., 74., 94., 80., 83., 70., 78., 85., 85., 56., 59., 56., 73., 33., 49., 50., 40., 22., 37., 26.,
6., 11., 7., 26., 0., 3., 0., 0., 0., 0., 0., 3., 9., 0., 31., 0., 11., 0., 8., 0., 9., 18.,9., 14., 0., 0., 6., 0.])

xdata = np.arange(0, len(ydata), 1)

plt.plot(xdata, ydata, 'bo', label='data')

popt, pcov = curve_fit(gauss, xdata, ydata, p0=[ydata.max(), ydata.mean(), ydata.std()])
plt.plot(xdata, gauss(xdata, *popt), 'r-', label='fit')

plt.show()

Solution 2:

You may not have really wanted to use ydata.mean() for the initial value of the Gaussian centroid, or ydata.std() for the initial value of the variance -- those are probably better guessed from xdata. I don't know if that is what caused the initial trouble.

You might find the lmfit library useful. This provides a way to turn your model function gauss into a Model class with a fit() method that uses named Parameters determined from your model function. Using it, your fit might look like:

import numpy as np
import matplotlib.pyplot as plt

from lmfit import Model

def gauss(x, a, x0, sigma):
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

ydata = np.array([2., 2., 11., 0., 5., 7., 18., 12., 19., 20., 36., 11., 21., 8., 13., 14., 8., 3., 21., 0., 24., 0., 12., 0., 8., 11., 18., 0., 9., 21., 17., 21., 28., 36., 51., 36., 47., 69., 78., 73., 52., 81., 96., 71., 92., 70., 84.,72., 88., 82., 106., 101., 88., 74., 94., 80., 83., 70., 78., 85., 85., 56., 59., 56., 73., 33., 49., 50., 40., 22., 37., 26., 6., 11., 7., 26., 0., 3., 0., 0., 0., 0., 0., 3., 9., 0., 31., 0., 11., 0., 8., 0., 9., 18.,9., 14., 0., 0., 6., 0.])

xdata = np.arange(0, len(ydata), 1)

# wrap your gauss function into a Model
gmodel = Model(gauss)
result = gmodel.fit(ydata, x=xdata, 
                    a=ydata.max(), x0=xdata.mean(), sigma=xdata.std())

print(result.fit_report())

plt.plot(xdata, ydata, 'bo', label='data')
plt.plot(xdata, result.best_fit, 'r-', label='fit')
plt.show()

There are several additional features. For example, you might want to see the confidence in the best fit, which would be (in master version, soon-to-be-released):

# add estimated band of uncertainty:
dely = result.eval_uncertainty(sigma=3)
plt.fill_between(xdata, result.best_fit-dely, result.best_fit+dely, color="#ABABAB")
plt.show()

to give: