试图用 scipy 拟合高斯函数

Question

我有一个衍射图，我需要将给定的峰拟合到高斯函数，我正在尝试使用

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

所以我做到了，但它在我的情况下不起作用

为了复现，可以使用这个数据集

检查第一列中的值是否接近另一个，我不知道是否与此有关，例如减法抵消，将e取负的小幂之类的；我找到了一些关于那个的东西，但我真的不明白

Answer 1

是的，你需要初始参数；不，你没有正确计算它们。 Y 的均值不是 随机变量的均值，Y 的标准差是not 随机变量的标准差。

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
mu0 = X[Y.argmax()]
sigma0 = 0.5
b0 = Y[:10].mean()
a0 = (Y.max() - 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, X.min(),                 0,       0),
        (np.inf, X.max(), X.max() - X.min(), Y.max()),
    ),
)
print(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()

你的数据也太吵了，没有足够的样本

m

没有任何意义。