python - 尝试用 scipy 拟合高斯

Question

我有一个衍射图，我需要将给定的峰值拟合到高斯函数，我尝试使用 scipy.optimize 中的 curve_fit 来实现这一点，但我出现一些错误。首先，由于我的数据，我在方程中添加了一条直线来拟合，所以我有

当我尝试通过以下方式拟合这个方程时

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 取负的小幂或其他什么；我发现了一些相关内容，但实际上我不明白

Answer 1

是的，你需要初始参数；不，你没有正确计算它们。 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
xmin, xmax = X.min(), X.max()  # left and right bounds
i_max = Y.argmax()             # index of highest value - for guess, assumed to be Gaussian peak
ymax = Y[i_max]     # height of guessed peak
mu0 = X[i_max]      # centre x position of guessed peak
b0 = Y[:20].mean()  # height of baseline guess

# https://en.wikipedia.org/wiki/Gaussian_function#Properties
# Index of first argument to be at least halfway up the estimated bell
i_half = np.argmax(Y >= (ymax + b0)/2)
# Guess sigma from the coordinates at i_half. This will work even if the point isn't at exactly
# half, and even if this point is a distant outlier the fit should still converge.
sigma0 = (mu0 - X[i_half]) / np.sqrt(
    2*np.log(
        (ymax - b0)/(Y[i_half] - b0)
    )
)

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.43599321e+02 6.86947536e+01 2.23407054e-01 3.02000000e+02]
Fit:   [2.28882618e+02 6.86885140e+01 2.25575284e-01 3.02398123e+02]

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