Python Spinmob curve_fit works but fitter does not

Question

I'm trying to fit data with a Gaussian.

The raw data itself displays a very obvious peak. When I attempt fitting using curve_fit, the fit identifies the peak but it does not have a curved top. I am trying to fit the data now with spinmob's fitter as well. However, this fitting just gives a straight line.

I've tried changing several parameters of the fitter, the Gaussian function definition, and the initial parameters for the fit but nothing seems to work.

Here is the code:

from scipy.optimize import curve_fit
from scipy import asarray as ar,exp
import spinmob as s

x = x30
y = ydata

def gaussian(x, A, mu, sig): # See http://mathworld.wolfram.com/GaussianFunction.html
    return A/(sig * np.sqrt(2*np.pi)) * np.exp(-np.power(x-mu, 2) / (2 * np.power(sig, 2)))

popt,pcov = curve_fit(gaussian,x,y,p0=[1,7.688,0.005])

FWHM = 2*np.sqrt(2*np.log(2))*popt[2]
print("FWHM: {}".format(FWHM))

plt.plot(x,y,'bo',label='data')
plt.plot(x,gaussian(x,*popt),'r+-',label='fit')
plt.legend()

fitter = s.data.fitter()
fitter.set(subtract_bg=True, plot_guess_zoom=True)
fitter.set_functions(f=gaussian, p='A=1,mu=8.688,sig=0.001')
fitter.set_data(x, y, eydata = 0.03)
fitter.fit()

The curve_fit returns this plot: Curve_fit plot

The spinmob fitter plot gives this: Spinmob Fitter Plot

As a test, if you use the values from spinmob as the p0 values for curve_fit and it works OK, then the problem should be the initial parameters being passed to curve_fit. would you please make this test? — James Phillips

M Newville M Newville · Accepted Answer · 2018-12-04T01:30:11

Assuming that spinmob actually uses scipy.curve_fit under the hood, I would guess (sorry) that the problem is that the initial values you give to it are so far off that it cannot possibly find a solution.

For sure, A=1 is not a very good guess for either scipy.curve_fit() or spinmob.fitter(). The peak is definitely negative, and you should be guessing a value more like -0.1 than +1. In fact you could probably assert that A must be < 0.

The initial value of 7.688 for mu that you give to curve_fit() is pretty good, and will allow a solution. I do not know whether it is a typo or not, but the initial value of 8.688 for mu that you give to spinmob.fitter() is very far off (that is, way outside the data range), and the fit will never be able to refine its way to the correct solution from there.

Initial values matter for curve-fitting and poor initial values can lead to bad results.

It might be viewed by some as a shameless plug, but allow me to encourage you to try lmfit (https://lmfit.github.io/lmfit-py/) (I am a lead author) for this kind of problem. Lmfit replaces the array of parameter values with named Parameter objects for better organization of fits. It also has a built-in Gaussian model (which also calculates FWHM, including an uncertainty). That is, with Lmfit, your script might look like:

import numpy as np
import matplotlib.pyplot as plt

from lmfit.models import GaussianModel
from lmfit.lineshapes import gaussian

# create fake data that looks like yours
xdata = 7.670 + np.arange(41)*0.0010
ydata = gaussian(xdata, amplitude=-0.196, center=7.6881, sigma=0.001)
ydata += np.random.normal(size=41, scale=10.0)

# create gaussian model
gmodel = GaussianModel()

# fit data, giving initial values for amplitude, center, and sigma
result = gmodel.fit(ydata, x=xdata, amplitude=-0.1, center=7.688, sigma=0.005)

# show results
print(result.fit_report())
plt.plot(xdata, ydata, 'bo', label='data')
plt.plot(xdata, result.best_fit, 'r+-', label='fit')
plt.legend()
plt.show()

This will print out a report like

[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 21
    # data points      = 41
    # variables        = 3
    chi-square         = 5114.87632
    reduced chi-square = 134.602009
    Akaike info crit   = 203.879794
    Bayesian info crit = 209.020510
[[Variables]]
    sigma:      9.7713e-04 +/- 1.5456e-04 (15.82%) (init = 0.005)
    center:     7.68822727 +/- 1.5484e-04 (0.00%) (init = 7.688)
    amplitude: -0.19273945 +/- 0.02643400 (13.71%) (init = -0.1)
    fwhm:       0.00230096 +/- 3.6396e-04 (15.82%) == '2.3548200*sigma'
    height:    -78.6917624 +/- 10.7894236 (13.71%) == '0.3989423*amplitude/max(1.e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
    C(sigma, amplitude) = -0.577

and produce a plot of data and best fit like

which should be close to what you are trying to do.

Python Spinmob curve_fit works but fitter does not

1 Answers