Search code examples
pythoncovariance

gaussian fit..error: OptimizeWarning: Covariance of the parameters could not be estimated


I am trying to fit my data with a gaussian function. I have set the parameters right and as far as I can evaluate, my code is also correct. But I am not getting the correct fit, and there is some error about the covariance of the parameters. Could anyone pls invest their time to review this code and tell me what I am missing?

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
old_settings = np.seterr(all='ignore')
from scipy import interp
from scipy import genfromtxt
import scipy
from math import pi, sqrt
old_settings = np.seterr(all='ignore')


data= genfromtxt('steel.txt')
x= data[:,0]
y= data[:,3]


def gaus(x, amp, mu, s):
    return amp *np.exp(-(x-mu)**2/(2*s**2))

amp,mu,s= 400, 0, -0.1
popt, pcov = curve_fit(gaus,x,y,p0=(400, 0, -0.1))
print(popt) 

p1 = popt[0] 
p2 = popt[1] 
p3 = popt[2] 



residuals = y - gaus(x, amp, mu, s)
fres = sum( (residuals**2)/gaus(x, amp, mu, s) ) 
print(fres)


curvey = gaus(x, p1,p2,p3)
plt.plot(x,y, 'b.', label='experimental data')
plt.plot(x, curvey, 'r.', label='gaussian fit')
plt.legend(loc='best')
plt.ylim(2700,4000)
plt.xlabel('velocity')
plt.ylabel('counts per seconds')
plt.legend()
plt.show()

My data is here: https://www.dropbox.com/s/7wn34goicl8wu0l/steel.txt?dl=0


Solution

  • Your fit function has a range from 0 to amp. This is not the range your data set has. Solution: Add an offset to your function:

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.optimize import curve_fit
    from scipy import genfromtxt
    
    data= genfromtxt('steel.txt')
    x= data[:,0]
    y= data[:,3]
    #added offset parameter
    def gaus(x, amp, mu, s, offset):
        return amp *np.exp(-(x-mu)**2/(2*s**2)) + offset
    
    popt, pcov = curve_fit(gaus,x,y,p0=(-4, 0, -0.1, 100))
    print(popt) 
    #better resolution for fit curve representation
    curvex = np.linspace(np.min(x), np.max(x), 1000)
    curvey = gaus(curvex, *popt)
    plt.plot(x,y, 'b.', label='experimental data')
    plt.plot(curvex, curvey, 'r', label='gaussian fit')
    plt.legend(loc='best')
    plt.ylim(2700,4000)
    plt.xlabel('velocity')
    plt.ylabel('counts per seconds')
    plt.legend()
    plt.show()
    

    Output enter image description here