试着用scipy拟合高斯分布

lf5gs5x2 于 2023-04-30 发布在其他

关注(0)|答案(1)|浏览(191)

我有一个衍射图，我需要将一个给定的峰拟合到一个高斯函数，我试图使用curve_fit，从scipy.optimize，但我得到了一些错误。首先，因为我的数据，我在方程中加了一条直线来拟合，所以我有

当我试着通过

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

def gauss(x, mu, sigma, m, b):
    return (1.0/(sigma*np.sqrt(2.0*np.pi)))*np.exp((-(x-mu)**2)/(2.0*sigma**2)) + m*x + b

#load the dataset
datadrx = np.genfromtxt('data.txt')

X = datadrx[:,0]
Y = datadrx[:,1]

mu0 = np.mean(Y)
sigma0 = np.std(Y)

popt, pcov = curve_fit(gauss, X, Y, p0 = [mu0,sigma0, 1.0, 1.0])
plt.title('Gaussian fit')
plt.plot(X, gauss(X,*popt))
plt.plot(X,Y)

我得到了这个警告

/home/eariasp/miniconda3/envs/herrcomp/lib/python3.10/site-packages/scipy/optimize/_minpack_py.py:906: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

和一条直线作为拟合，看到了吗

当我移除mx+b部分时，我得到了“更好”的匹配（我的意思是更好，因为至少它不仅仅是一条直线），但它仍然远离我想要的匹配

所以，我的问题是，什么是错的？我发现定义猜测参数p0是个好主意，所以我就这么做了，但它在我的例子中不起作用
为了重现，可以使用this数据集
检查第一列的值是否接近，我不知道这是否有关系，比如减法抵消，取e为负的小幂或其他什么;我发现了一些关于那个的东西，但我真的没有得到它

scipy

来源：https://stackoverflow.com/questions/76137714/trying-to-fit-a-gaussian-with-scipy

1条答案

按热度按时间

jdgnovmf1#

是的，你需要初始参数;不，你计算得不对。Y的均值不是随机变量的均值，Y的标准差不是随机变量的标准差。
约束条件和猜测越好，拟合收敛的速度就越快、越可靠。值得庆幸的是，对于高斯模型，这样的约束和猜测是相当简单的计算-事实上，猜测是 * 如此 * 好，对于一些（不是所有）应用程序，甚至不需要拟合。

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

def gauss(x: np.ndarray, a: float, mu: float, sigma: float, b: float) -> np.ndarray:
    return (
        a/sigma/np.sqrt(2*np.pi)
    )*np.exp(
        -0.5 * ((x-mu)/sigma)**2
    ) + b

X, Y = np.loadtxt('data.txt', delimiter=' ').T
yimax = Y.argmax()
ymax = Y[yimax]
xmin, xmax = X.min(), X.max()

mu0 = X[yimax]
b0 = Y[:10].mean()

# https://en.wikipedia.org/wiki/Gaussian_function#Properties
half_width = mu0 - X[np.argmax(Y >= (ymax + b0)/2)]
sigma0 = half_width / np.sqrt(2*np.log(2))

a0 = (ymax - b0) * sigma0 * np.sqrt(2*np.pi)
p0 = a0, mu0, sigma0, b0

popt, _ = curve_fit(
    f=gauss, xdata=X, ydata=Y, p0=p0,
    bounds=(
        (     1, xmin,           0,    0),
        (np.inf, xmax, xmax - xmin, ymax),
    ),
)
print('Guess:', np.array(p0))
print('Fit:  ', popt)

fig, ax = plt.subplots()
ax.set_title('Gaussian fit')
ax.scatter(X, Y, marker='+', label='experiment', color='orange')
ax.plot(X, gauss(X, *p0), label='guess', color='lightgrey')
ax.plot(X, gauss(X, *popt), label='fit')
ax.legend()
plt.show()

Guess: [2.45543591e+02 6.86947536e+01 2.23036703e-01 2.97800000e+02]
Fit:   [2.28882519e+02 6.86885141e+01 2.25575115e-01 3.02398139e+02]

你的数据也太嘈杂，没有足够的样本让m有任何意义。

赞(0）回复(0）举报 2023-04-30

我来回答

试着用scipy拟合高斯分布

1条答案

相关问题

热门标签

最新问答