来自 scipy.optimize 的曲线拟合问题答案

【问题标题】：Problems with curve_fit from scipy.optimze来自 scipy.optimize 的曲线拟合问题
【发布时间】：2016-08-20 00:41:27
【问题描述】：

我知道有一些类似的问题，但由于它们都没有让我更进一步，我决定问我自己的一个。很抱歉，如果我的问题的答案已经在某个地方，但我真的找不到。

我尝试使用 curve_fit 将 f(x) = a*x**b 拟合到相当线性的数据。它编译正确，但结果如下所示：

问题是，我真的不知道自己在做什么，但另一方面，试穿更像是一门艺术而不是科学，而且至少有一位将军 bug with scipy.optimize。

我的数据如下所示：

x 值：

[16.8, 2.97, 0.157, 0.0394, 14.000000000000002, 8.03, 0.378, 0.192, 0.0428, 0.029799999999999997, 0.000781, 0.0007890000000000001]

y 值：

[14561.766666666666, 7154.7950000000001, 661.53750000000002, 104.51446666666668, 40307.949999999997, 15993.933333333332, 1798.1166666666666, 1015.0476666666667, 194.93800000000002, 136.82833333333332, 9.9531566666666684, 12.073133333333333]

这是我的代码（在that question 的最后一个答案中使用了一个非常好的示例）：

def func(x,p0,p1): # HERE WE DEFINE A FUNCTION THAT WE THINK WILL FOLLOW THE DATA DISTRIBUTION
    return p0*(x**p1)

# Here you give the initial parameters for p0 which Python then iterates over to find the best fit
popt, pcov = curve_fit(func,xvalues,yvalues, p0=(1.0,1.0))#p0=(3107,0.944)) #THESE PARAMETERS ARE USER DEFINED

print(popt) # This contains your two best fit parameters

# Performing sum of squares
p0 = popt[0]
p1 = popt[1]
residuals = yvalues - func(xvalues,p0,p1)
fres = sum(residuals**2)

print 'chi-square'
print(fres) #THIS IS YOUR CHI-SQUARE VALUE!

xaxis = np.linspace(5e-4,20) # we can plot with xdata, but fit will not look good 
curve_y = func(xaxis,p0,p1)

起始值来自与 gnuplot 的匹配，这是合理的，但我需要交叉检查。

这是打印输出（首先拟合 p0、p1，然后是卡方）：

[  4.67885857e+03   6.24149549e-01]
chi-square
424707043.407

我想这是一个很难的问题，因此提前非常感谢！

【问题讨论】：

标签： python matplotlib curve-fitting

【解决方案1】：

拟合curve_fit时优化(data - model)^2 / (error)^2之和

如果您不传递错误（就像您在此处所做的那样）curve_fit 假定所有点的错误均为 1。

在这种情况下，由于您的数据跨越多个数量级，具有最大 y 值的点支配目标函数，并导致 curve_fit 尝试以牺牲其他点为代价来拟合它们。

解决此问题的最佳方法是将您的yvalues 中的错误包含在合适的位置（看起来您这样做，因为您制作的情节中有错误栏！）。您可以通过将它们作为curve_fit 的sigma 参数传递来做到这一点。

【讨论】：

非常感谢，使用sigma 让我看起来很合身。我以前没有使用过它，因为我认为它 a) 不重要并且 b) 用于 xvalues 中的错误。（我不明白为什么xvalues 中没有错误参数。）此外，我认为估计错误为 1 是设计缺陷，更喜欢0.1*value 之类的东西。你会同意吗？你认为值得提交错误报告吗？问候。
不是真的 - 在代码中对我来说最不意外的路径是，如果没有给出错误，那么它应该假设错误是恒定的。这通常是你想要的。 documentation is also pretty clear here
我不确定。我一直认为，它必须以一种对点的重要性没有影响的方式，因为先验它们之间没有区别。你知道，0是否可能？

【解决方案2】：

我会重新考虑实验部分。有两个数据点值得怀疑：

你给我们看的图片看起来不错，因为你拿了日志：

您可以对 log(x) 和 log(y) 进行线性拟合。通过这种方式，您可能会限制最大残差的影响。另一种方法是稳健回归（来自 sklearn 的 RANSAC 或来自 scipy 的 minimum_squares）。

不过，您应该收集更多数据点或重复测量。

【讨论】：

好吧，其他数据并不是一个真正的选择，但如果不合理，我可以忽略点，因为我所有的点都是单一的错误。但是您对计算的想法非常有帮助，我会记住它们。谢谢。