英文:
How to fix gaussian fit not behaving like expected?
问题
我有一组数据，显示辐射的吸收与速度的函数关系。数据显示一个非常明显的低谷，或者如果我们绘制数据的倒数，即吸收度，我们会得到一个明显的峰值。我没有理由不相信这个峰值是一个高斯分布，想要进行拟合以获取这个峰值的方差。因此，我尝试使用scipy.optimize.curve_fit来实现这一目标，无论初始猜测如何，得到的拟合结果都相差很远。
我附上了代码和结果图：
import matplotlib.pyplot as pltfrom scipy.optimize import curve_fitimport numpy as np# 以下是数据的读取和处理部分，未翻译# ...def func(x, mu, sigma):return (1 / (np.sqrt(np.pi * 2) * sigma)) * np.exp(-0.5 * ((x - mu) / sigma) ** 2)data = np.array(cts)cts_norm = (data - data.min()) / (data.max() - data.min())cts_inv = 1 - cts_normfit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2])print(fit)plt.plot(vel, cts_inv, 'bo')plt.plot(vel, func(vel, fit[0], fit[1]), 'r')
请注意，以上是代码的翻译，没有包含问题的回答或其他内容。
英文:
I have a set of data showing radition not being absorbed as a function of velocity. The data shows a very clear dip or if we plot the inverse of the data the absorbtion, we get a clear peak instead. I have no reason not to belive this peak to be a gaussian and would like to make a fit to get the variance of this peak. So I've tried to use scipy.optimize.curve_fit, to achieve this. Both using scipy.stats.norm.pdf and a self written version of the function. No matter initial guesses. The resulting fit is way of.
I attached the code and a picture of the resulting graph.
What am I doing wrong? Are there any general tricks for these kind of tasks?
import matplotlib.pyplot as pltfrom scipy.optimize import curve_fitimport numpy as npcts = []vel = []file = open("Heisenberg/Mössbauer/Final.lst", "r")linesArr = file.readlines()for i in range(210, 260):lineList1 = linesArr[i].split()cts.append(int(lineList1[1]))chn = (int(lineList1[0]))tempVel = -0.04 * chn + 9.3vel.append(tempVel)def func (x, mu,sigma):return (1 / (np.sqrt(np.pi * 2) *sigma)) * np.exp(-0.5*((x-mu)/sigma)**2)data = np.array(cts)cts_norm = (data - data.min())/ (data.max() - data.min())cts_inv = 1 - cts_normfit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2])print(fit)plt.plot(vel, cts_inv, 'bo')plt.plot(vel, func(vel, fit[0],fit[1]),"r")
答案1
得分: 2
问题在于您试图拟合一个正态分布到不是概率分布的数据上！概率分布的积分等于1，但对于您的数据来说并非如此，它可以具有任何幅度。很难将您的数据归一化以满足这一要求。相反，您可以简单地添加一个新的参数来控制正态分布的“幅度”，如下所示：
import matplotlib.pyplot as pltfrom scipy.optimize import curve_fitimport numpy as npcts = [0, 0, 0, 0, -1, -2, -5, -10, -5, -2, -1, 0, 0, 0, 0]vel = np.linspace(-0.75, 1.25, 15)def func(x, mu, sigma, a):return a * np.exp(-0.5 * ((x - mu) / sigma) ** 2)data = np.array(cts)cts_norm = (data - data.min()) / (data.max() - data.min())cts_inv = 1 - cts_normfit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1])print(fit)plt.plot(vel, cts_inv, 'bo')plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r")plt.show()
（我使用了一些虚拟数据，因为我无法访问您的文件，但这并不重要）
英文:
The issue is that you are trying to fit a normal distribution with data that is not a probability distribution! Probability distributions have an integral equal to 1, but that is not the case for your data, which can have any amplitude. It would be hard to normalize your data to satisfy this. Instead, you can simply add a new parameter which controls the "amplitude" of the normal distribution, like below:
import matplotlib.pyplot as pltfrom scipy.optimize import curve_fitimport numpy as npcts = [0, 0, 0, 0, -1, -2, -5, -10, -5, -2, -1, 0, 0, 0, 0]vel = np.linspace(-0.75, 1.25, 15)def func(x, mu, sigma, a):return a * np.exp(-0.5 * ((x - mu) / sigma) ** 2) # << heredata = np.array(cts)cts_norm = (data - data.min()) / (data.max() - data.min())cts_inv = 1 - cts_normfit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1]) # << hereprint(fit)plt.plot(vel, cts_inv, 'bo')plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r") # << and hereplt.show()
(I used some dummy data as I don't have access to your file, but it doesn't really matter)
答案2
得分: 1
我会为你的模型添加一些灵活性，具体如下。我通过拍摄图像并使用 [此][1] 免费网络服务来检索了你的数据。```pythonimport matplotlib.pyplot as pltimport numpy as npfrom scipy.optimize import curve_fitfrom scipy.stats import normdata = np.loadtxt(r"C:\Users\Cristiano\Desktop\data.txt", delimiter=",")x = data[:, 0]y = data[:, 1]def f(x, a, b, mu, sigma):return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))popt, pcov = curve_fit(f, x, y)mean, std = norm.fit(data)plt.scatter(x, y)xx = np.linspace(-0.75, 1.25, 1000)plt.plot(xx, f(xx, *popt))plt.show()
<details><summary>英文:</summary>I would add a little more flexibility to your model as follows. I retrieved your data by taking a screenshot of the image and using [this][1] free web service.```pythonimport matplotlib.pyplot as pltimport numpy as npfrom scipy.optimize import curve_fitfrom scipy.stats import normdata = np.loadtxt(r"C:\Users\Cristiano\Desktop\data.txt", delimiter=",")x = data[:, 0]y = data[:, 1]def f(x, a, b, mu, sigma):return a + b * np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))popt, pcov = curve_fit(f, x, y)mean, std = norm.fit(data)plt.scatter(x, y)xx = np.linspace(-0.75, 1.25, 1000)plt.plot(xx, f(xx, *popt))plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论