scipy curve_fit does not find the best fit

I'd like to find out whether there is a peak in my time series or not. To do so, I'm trying to fit the data points to a gauss curve. It works well for thousands of my samples:

but few are not fitter properly despite there is an obvious peak:

(see the very low peak with the highest point around 0.03)

Here's the code:

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

param, covariance = curve_fit(gauss, x, y)

I've noticed that the magnitude of the y values plays a role in the fitting so I've rescaled the data into <0, 100> interval. This helped but did not solve all the cases. Is there anything else I can do to improve the fitting? Different initialization? Smaller optimization step?

Here are some facts about the data:

• Every sample has 3-20 data points.
• The peak (if there is any) must have its highest point inside the span.
• x axis is from 0 to 20
• y axis is from 0 to 100

I've browse through other similar questions at stackoverflow but did not find a solution to my problem.

If anybody knows a better solution to determine whether there is a peak in the time series or not, I'd appreciate to hear it in the comments. Anyway, I'd like to know why some curves are not fitted properly.

The nonlinear regression involves iterative calculus starting from initial "guessed" values of the parameters.

Probably the trouble appears when the "guessed" values are too far from the unknown correct values.

In order to avoid the difficulty one can use a non-iterative method which doesn't requires initial values.

This non-conventional method is shown below :

Numerical example :

The general principle of this method is explained in : https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales . Instead of fitting the model equation one fit an integral equation to which the model equation is solution. This requires preliminary numerical integrations (for Sk and Tk). If they are not enough points the numerical integrations are not accurate enough and the final fitting is bad. That is why I didn't try to fit the first OP's example : A data of 6 points is obviously not enough. The non-iterative method of regression with integral equation is not recommended in this case.

Note : The coordinates of the points were measured on the second OP's graph which is not accurate measurement.

I suggest to the OP to try the above resuts as initial values of the parameters to start the iterative process of his software and see if the trouble of calculus disappears. Taking care of the different symbols used : his "a"= my "b" ; his "x0"= my "c" ; his "sigma"= my "1/sqrt(2p)".

The curve_fit method provides an option to define the initial function parameters. Choosing the initial height and position of the bell curve to the y and x coordinates of the "highest" point in the series works very well on my data.

Additionally, curve_fit also allows to set bounds of the function parameters to assure the height of the bell curve does not go below (or above) certain level.

Following code works for all of my data samples:

init = [max(y), np.argmax(y), 1.5]
bounds = ([50, len(x) * 0.05, -20], [150, (len(x) - 1) * 0.95, len(x) / 2])
param, covariance = curve_fit(gauss, x, y, p0=init, bounds=bounds)

If the parameters are outside the bounds, I catch the exception (ValueError) and assume there is no peak in the time series.

The Polynomial Feature will return a curve set for each degree of freedom. I demonstrated a 3rd order polynomial with two of its curve fits.

from sklearn.preprocessing import PolynomialFeatures

def gaussian(x, mu, sig):
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

n=500
x = np.linspace(0, 1, n)
mu=x.mean()
sig=x.std()

y=gaussian(x,mu,sig)

poly = PolynomialFeatures(degree=3)

curve_y = poly.fit_transform(np.array(y).reshape(-1,1))
print(curve_y)

#three degrees of free
degree_of_freedom=1
plt.plot(x, curve_y[:,degree_of_freedom], '.',c="red")
degree_of_freedom=2
plt.plot(x, curve_y[:,degree_of_freedom], '.',c="red")
plt.plot(x, y, '.',alpha=0.2, c="blue")
plt.show()

print("you can see degree of 1 is the best")