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
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()