Python - 不合适答案

【问题标题】：Python - wrrong fitPython - 不合适
【发布时间】：2020-01-14 13:14:21
【问题描述】：

我正在尝试重现已知的拟合结果（在期刊论文中报告）：将幂律模型应用于数据。从下面的图 A 中可以看出，我能够通过使用已知的最佳拟合参数来重现结果。

但是，我无法自己重新推导最适合的参数。

案例-A 返回，

OptimizeWarning: Covariance of the parameters could not be estimated（如果我省略了几个初始数据点，拟合会返回一些不错的结果，但仍与已知的最佳拟合结果不同）。 编辑： 现在，我这次只是发现了新的附加错误消息..： (1)RuntimeWarning: overflow encountered in power (2)RuntimeWarning: invalid value encountered in power

case-B（初始猜测更接近最佳拟合参数）返回，

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 5000. 如果我将 maxfev 设置为更高以考虑此错误消息，则拟合有效但返回错误结果（与最佳拟合结果相比非常错误）。

【问题讨论】：

请提供更完整的脚本以适合您，并附上结果。恐怕这是一个非常常见的问题，但是：由于您将数据绘制在对数图上，您是否也适合 log(y) 和 log(x)？由于您的 y 数据变化 5 或 6 个数量级，如果您不适合日志空间，则只有具有最高 y 值的 3 或 4 个数据点才重要。
当我直接使用“最佳拟合”参数在对数刻度上绘图时，绘图与发布的图像相似。然而，在普通最小二乘 (OLS) 回归中，回归误差的总和为零 - 但这里的情况并非如此。我建议重复这个测试，使用“最佳拟合”值计算回归误差的总和。如果您的测试中也出现这种情况，那么可能是作者没有使用 OLS。我还尝试使用“最佳拟合”参数作为曲线拟合的初始参数估计值，但奇怪的是这并没有收敛——我建议也进行该测试。
@MNewville，感谢您的评论。我刚刚用我使用的详细脚本更新了帖子。之前尝试过删除几个初始数据点，但似乎不是原因。
@James，感谢您的评论。我可以问你如何执行你建议的测试吗？如果不是 OLS，那么正确匹配的替代方法是什么？
要确定回归误差，请将“最佳拟合”参数直接传递给bendp() 函数，并将您的数据作为“x”以找到预测的“最佳拟合”值。在这种情况下，“最佳拟合”回归的预测误差将是“预测的 - y”，如果使用 OLS，这些误差的总和应为零。如果总和不为零，我的建议是联系论文的作者，询问他们是否使用 OLS 对您拥有的数据进行回归，如果没有，他们是如何进行回归的。

标签： python curve-fitting covariance power-law loglog

【解决方案1】：

正如我评论的那样：

由于您在对数图上绘制数据，您是否也适合日志（y）和日志（x）？由于您的 y 数据变化 5 或 6 个订单幅度，如果您不适合对数空间，则只有那 3 或 4 个数据 y 值最高的点很重要。

显然这暗示有点太微妙了。所以我会更直接：如果您在日志空间中绘图，则适合日志空间。

而且，您的模型很容易从 negative**fraction 和 NaN 生成复数，这无疑会导致您的所有拟合出现问题。总是打印出最佳拟合参数和协方差矩阵。

因此，您可能需要对参数施加限制（当然，我不知道您的模型是否正确，或者实际上是您认为“正确答案”使用的）。也许从这样的事情开始：

import matplotlib.pyplot as plt
from lmfit import Model
import numpy as np

# the data can be found at the end of this post.
xx, yy = np.genfromtxt('the_data', unpack=True)

# the BPL function
def bendp(x, A, x_bend, allo, alhi, c):
    numer = A * x**-allo
    denom = 1 + (x/x_bend)**(alhi-allo)
    out = (numer/denom) + c
    return  np.log(out)       ## <- TAKE THE LOG


mod = Model(bendp)
para = mod.make_params(A = 0.01, x_bend=1e-2, allo=1., alhi=2., c=1e-2)

# Note: your model is very sensitive # make sure A, x_bend, alhi, and c cannot be negative
para['A'].min = 0
para['x_bend'].min = 0
para['alhi'].min = 0
para['alhi'].max = 3
para['c'].min = 0


result = mod.fit(np.log(yy), para, x=xx) ## <- FIT THE LOG

print(result.fit_report())

plt.loglog(xx, yy, linewidth=0, marker='o', markersize=6)
plt.loglog(xx, np.exp(result.best_fit), 'r', linewidth=5)
plt.show()

希望对您有所帮助...

【讨论】：