Fitting with a gaussian

Question

I have some problems when trying to fit data from a text file with a gaussian. This is my code, where cal1_p1 is an array containing 54 values.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np  
from scipy.optimize import curve_fit

cal1=np.loadtxt("C:/Users/Luca/Desktop/G3/X_rays/cal1_5min_Am.txt")
cal1_p1=[0 for a in range(854,908)]
for i in range(0,54):
      cal1_p1[i]=cal1[i+854]   
# cal1_p1 takes the following values:

[5.0,6.0,5.0,11.0,4.0,9.0,14.0,13.0,13.0,14.0,12.0,13.0,16.0,20.0,15.0,23.0,23.0,33.0,43.0,46.0,41.0,40.0,49.0,57.0,62.0,61.0,53.0,65.0,64.0,42.0,72.0,55.0,47.0,43.0,38.0,46.0,37.0,39.0,27.0,18.0,20.0,20.0,18.0,10.0,11.0,8.0,10.0,6.0,8.0,8.0,6.0,10.0,6.0,4.0]

x=np.arange(854,908)

def gauss(x,sigma,m):
    return np.exp(-(x-m)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))
from scipy.optimize import curve_fit
popt,pcov=curve_fit(gauss,x,cal1_p1,p0=[10,880])

plt.xlabel("Channel")
plt.ylabel("Counts")
axes=plt.gca()
axes.set_xlim([854,907])
axes.set_ylim([0,75])
plt.plot(x,cal1_p1,"k")     
plt.plot(x,gauss(x,*popt),'b', label='fit')

The problem is that the resulting gaussian is really squeezed, namely it has a very low variance. Even if I try to modify the initial value p_0 the result doesn't change. What could be the problem? Thanks for any help you can provide!

What is `cal1_p1 `? Share the working code so that your result can be reproduced. — Sheldore, Aug 16 '18 at 08:27
Possibly a duplicate. Does the answer in https://stackoverflow.com/questions/19206332/gaussian-fit-for-python help? — Ghasem Naddaf, Aug 16 '18 at 08:30
@GhasemNaddaf If I try to do that, it tells me : "Covariance of the parameters could not be estimated". This is due to the fact that the resulting values of the mean and of sigma are really higher than they should be. — Luca, Aug 16 '18 at 08:49

score 1 · Accepted Answer · answered Aug 16 '18 at 08:54

The problem is that the Gaussian is normalised, while your data are not. You need to fit an amplitude as well. That is easy to fix, by adding an extra parameter a to your function:

x = np.arange(854, 908)

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

popt, pcov = curve_fit(gauss, x, cal1_p1, p0=[10, 880, 1])
print(popt)
plt.xlabel("Channel")
plt.ylabel("Counts")
axes=plt.gca()
axes.set_xlim([854, 907])
axes.set_ylim([0, 75])
plt.plot(x, cal1_p1, "k")  
plt.plot(x, gauss(x,*popt), 'b', label='fit')

While I've given 1 as starting parameter for a, you'll find that the fitted values are actually:

[   9.55438603  880.88681556 1398.66618699]

but the amplitude value here can probably be ignored, since I assume you'd only be interested in the relative strength, which can be measured in counts.

Fitting with a gaussian

1 Answers1