python - 为什么 SciPy 的curve_fit 找不到协方差/给我这个高阶高斯函数的感觉参数?
问题描述
这是一些最小的代码:
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] 供参考),但没有帮助。
解决方案
我相信您遇到了这个问题,因为curve_fit
也在测试 的非整数值n
,在这种情况下,您的函数gauss
在x<x_o
.
我相信通过暴力破解每个整数并为每个整数n
找到最佳参数会更容易。例如,如果您考虑最多 50 个,则 的选项很少,而暴力破解实际上是一种不错的方法。查看您的数据,最好只考虑2 的倍数(除非您的其他数据看起来可能是奇数)。sig,c,x_o
n
n = 0,1,2,...
n
n
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]
。
这是结果图:(蓝线是数据,橙线是高斯拟合)
推荐阅读
- javascript - 如何在 vue cli 项目中集成 trading-vue-js?
- javascript - Reactjs - 如何根据成功的 axios POST 创建布尔值?
- asp.net - 从另一个组件更新角度数据源会使角度数据表分页不起作用
- python - Matplotlib如何为直方图的子图添加全局图例
- python - How do you add an extra value to a dictionary comprehension with enumerate?
- r - 使用 dplyr R 删除在超过“n”个时间点具有零值的 ID
- node.js - Puppeteer 启动多个具有唯一数据的实例?
- html - 将网格模板设置为相同的宽度和高度
- python - 分类和数值数据的混合 pyplot
- google-sheets - 谷歌表格中的单元格验证日期