scipy 如何修复高斯拟合不像预期的行为?

3phpmpom  于 2023-01-13  发布在  其他
关注(0)|答案(2)|浏览(238)

我有一组数据,显示辐射没有被吸收为速度的函数。数据显示了一个非常明显的下降,或者如果我们绘制吸收数据的倒数,我们会得到一个清晰的峰。我没有理由不相信这个峰是高斯峰,并希望拟合得到这个峰的方差。所以我尝试使用scipy.optimize.curve_fit,来实现这个。使用scipy.stats.norm.pdf和一个自己编写的函数版本。不管最初的猜测。结果拟合是的方式。我附上了代码和结果图的图片。我做错了什么?有什么一般技巧,为这类任务?

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
cts = []
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.3
    vel.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_norm
fit, 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")

h22fl7wq

h22fl7wq1#

问题在于,您试图用非概率分布的数据拟合正态分布!概率分布的积分等于1,但您的数据却不是这样,它可以具有任何幅度。要使数据满足这一要求,将数据归一化是很困难的。相反,您可以简单地添加一个新参数来控制正态分布的“幅度”,如下所示:

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

cts = [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)  #      << here

data = np.array(cts)
cts_norm = (data - data.min()) / (data.max() - data.min())
cts_inv = 1 - cts_norm
fit, error = curve_fit(func, vel, cts_inv, p0=[0.2, 0.2, 1]) # << here
print(fit)

plt.plot(vel, cts_inv, 'bo')
plt.plot(vel, func(vel, fit[0], fit[1], fit[2]), "r")  #      << and here
plt.show()

(我使用了一些虚拟数据,因为我无法访问您的文件,但这并不重要)

i7uaboj4

i7uaboj42#

我想给你的模型增加一点灵活性,如下所示:我通过截取图像并使用this免费web服务来检索你的数据。

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import norm

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

相关问题