CODEWITHSUNDEEP
pythonscipycurve-fittinggaussian
James
James

Reputation: 67

scipy.optimize.curve_fit can't fit an upside down gaussian

I'm trying to fit a line to an upside down gaussian distribution using scipy.optimize.curve_fit. It works perfectly to fit a traditional gaussian, but wont fit a gaussian with the sign flipped, and instead will always output a straight line.

I've also tried constraining my gaussian function so that the variable 'a' is always negative but this does not solve the problem. Specifying -max(y) instead of max(y) also does not seem to help.

import scipy
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import asarray as ar,exp

def fitdata(x, y):
    mean = sum(x * y) / sum(y)
    sigma = np.sqrt(sum(y * (x - mean) ** 2) / sum(y))

    def guassianfunc(xVar, a, b, c):
        return a * exp(-(xVar - b) ** 2 / (2 * c ** 2))

    popt, _ = scipy.optimize.curve_fit(guassianfunc, x, y, p0=[max(y), mean,         sigma])
    return guassianfunc(np.arange(1, 6, 1), *popt)


x = np.array((1,2,3,4,5))
f, ((ax1, ax2)) = plt.subplots(2, sharex='col', sharey='row')
y = np.array((1, 2, 3, 2, 1))
ax1.plot(x, y, color='black')
ax1.plot(x, fitdata(x, y), linewidth=2, label='Fit')
y = np.array((3, 2, 1, 2, 3))
ax2.plot(x, y, color='black')
ax2.plot(x, fitdata(x, y), linewidth=2, label='Fit')
plt.legend()

output

Upvotes: 3

Views: 1518

Answers (1)

Dan
Dan

Reputation: 45752

I think it's because your second line would require an offset parameter maybe? i.e. your guassianfunc can make an upside down bell-curve provided it's below the x-axis. For example:

y = np.array((3, 2, 1, 2, 3))-4
ax2.plot(x, y, color='black')
ax2.plot(x, fitdata(x, y), linewidth=2, label='Fit')

Try adding in a 4th param like this

def fitdata(x, y):
    mean = sum(x * y) / sum(y)
    sigma = np.sqrt(sum(y * (x - mean) ** 2) / sum(y))

    def guassianfunc(xVar, a, b, c, d):
        return a * exp(-(xVar - b) ** 2 / (2 * c ** 2)) + d

    popt, _ = scipy.optimize.curve_fit(guassianfunc, x, y, p0=[max(y), mean,         sigma,0])
    return guassianfunc(x, *popt)

This also reduces your error significantly.

Upvotes: 3

Related Questions