为什么 SciPy 的 curve_fit 找不到协方差/给我这个高阶高斯函数的感觉参数？答案

【问题标题】：Why can't SciPy's curve_fit find the covariance/give me sensical parameters for this higher order gaussian function?为什么 SciPy 的 curve_fit 找不到协方差/给我这个高阶高斯函数的感觉参数？
【发布时间】：2021-08-25 02:39:10
【问题描述】：

这是一些最小的代码：

from scipy.optimize import curve_fit
xdata = [16.530468600170202, 16.86156794563677, 17.19266729110334, 17.523766636569913,
 17.854865982036483, 18.18596532750305, 18.51706467296962, 18.848164018436194, 19.179263363902763,
 19.510362709369332]
ydata = [394, 1121, 1173, 1196, 1140, 1196, 1158, 1160, 1046, 416]
#function i'm trying to fit the data to (higher order gaussian function)
def gauss(x, sig, n, c, x_o):
    return c*np.exp(-(x-x_o)**n/(2*sig**n))
popt = curve_fit(gauss, xdata, ydata)
#attempt at printing out parameters
print(popt)

当我尝试执行代码时，我收到以下错误消息：

ExCurveFit:9: RuntimeWarning: invalid value encountered in power
  return c*np.exp(-(x-x_o)**n/(2*sig**n))
C:\Users\dsneh\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.8_qbz5n2kfra8p0\LocalCache\local-packages\Python38\site-packages\scipy\optimize\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',
(array([nan, nan, nan, nan]), array([[inf, inf, inf, inf],
       [inf, inf, inf, inf],
       [inf, inf, inf, inf],
       [inf, inf, inf, inf]]))

我已经看到我应该忽略第一个，但也许这是一个问题。第二个更令人担忧，我显然希望参数值更有意义。我尝试添加参数作为猜测（[2,3,1000,17] 供参考），但没有帮助。

【问题讨论】：

看来优化过程找不到合适的最优值。您是否尝试过使用不同的methods 例如trf 并使用bounds 限制参数的值？
我没有。问题是对于某些 n 值的输出很复杂，所以如果可能的话，我想将其限制为 0.4 的倍数。这可以通过暴力破解来实现，从初步的谷歌看来，trf 似乎可以处理“稀疏”问题。这是否意味着我可以调整它来解决问题而无需暴力破解？
您可以使用bounds 属性将x_o 限制为始终小于最低数据点x。这将防止您的输出变得复杂。但是，bounds 不能用于将 n 限制为仅限整数。问题是您的数据是否可以用这样的x_o 来解释。在不知道您的实际数据的情况下，我猜x_o 肯定不应该受到限制。我认为使用蛮力n 方法你的成功机会应该更高。
这是有道理的。谢谢！

标签： python numpy scipy curve-fitting

【解决方案1】：

我相信您遇到了这个问题，因为curve_fit 也在测试n 的非整数值，在这种情况下，您的函数gauss 在x<x_o 时返回复数。

我相信暴力破解每个整数 n 并为每个 n 找到最佳参数 sig,c,x_o 会更容易。例如，如果你认为n = 0,1,2,... 可能高达 50，那么n 的选项就很少，而暴力破解实际上是一种不错的方法。查看您的数据，最好仅将 n 视为 2 的倍数（除非您的其他数据看起来像 n 可能是奇数）。

此外，您可能应该为其他参数引入一些界限，例如拥有sig>0 会很好，并且您可以安全地设置c>0（也就是说，如果您的所有数据看起来都像您的最小数据包含在您的问题中）。

这是我所做的：

p0 = [1,1200,18] # initial guess for sig,c,x_o
bounds = ([0.01,100,10],[1000,2500,200]) # (lower,upper) bounds for sig,c,x_o
min_err = np.inf # error to minimize to find the best value for n

for n in range(0,41,2): # n = 0,2,4,6,...,40

    def gauss(x, sig, c, x_o, n=n):
        return c*np.exp(-(x-x_o)**n/(2*sig**n))
    
    popt,pcov = curve_fit(gauss, xdata, ydata, p0=p0, bounds=bounds)

    err = np.sqrt(np.diag(pcov)).sum()
    if err < min_err:
        min_err, best_n, best_popt, best_pcov = err, n, popt, pcov

print(min_error, best_n, best_popt, best_pcov, sep='\n')

import matplotlib.pyplot as plt
plt.plot(xdata,ydata)
plt.plot(xdata, gauss(xdata, *best_popt, n=best_n))
plt.show()

我得到best_n = 10 和best_popt = [1.38, 1173.52, 18.02]。

这是结果图：（蓝线是数据，橙线是高斯拟合）

【讨论】：

我明白了，谢谢！我不知何故错过了输出有时很复杂的事实。