时间戳

首页 > 解决方案 > 用高斯拟合

python - 用高斯拟合

问题描述

尝试使用高斯拟合来自文本文件的数据时遇到一些问题。这是我的代码，其中 cal1_p1 是一个包含 54 个值的数组。

%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')

问题是得到的高斯真的被压缩了，即它的方差非常低。即使我尝试修改初始值 p_0，结果也不会改变。可能是什么问题呢？感谢您的任何帮助，您可以提供！

标签： pythonnumpy

解决方案

问题是高斯是标准化的，而您的数据没有。您还需要拟合振幅。通过向函数添加一个额外的参数，这很容易解决a：

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

虽然我给出了 1 作为 的起始参数a，但您会发现拟合值实际上是：

[   9.55438603  880.88681556 1398.66618699]

但是这里的幅度值可能会被忽略，因为我假设你只对相对强度感兴趣，它可以用计数来测量。

推荐阅读