在numpy和scipy中进行指数衰减曲线拟合

Question

我在拟合数据曲线时遇到了一些困难，但是却不知道出了什么问题。

在以前的工作中，我使用 numpy.linalg.lstsq 来拟合指数函数，而使用 scipy.optimize.curve_fit 来拟合Sigmoid函数。这次我希望创建一个脚本，让我可以指定各种函数，确定参数并测试其与数据的拟合情况。在这个过程中我注意到 Scipy 的 leastsq 和 Numpy 的 lstsq 对于相同的数据和相同的函数提供了不同的答案。这个函数很简单，是 y = e^(l*x) ，且在 x=0 时有 y=1 的限制。

Excel趋势线与Numpy的 lstsq 结果一致，但是由于Scipy的 leastsq 可以处理任何函数，因此最好弄清楚问题所在。

import scipy.optimize as optimize
import numpy as np
import matplotlib.pyplot as plt

## Sampled data
x = np.array([0, 14, 37, 975, 2013, 2095, 2147])
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,     0.001485394,     0.000495131])

# function
fp = lambda p, x: np.exp(p*x)

# error function
e = lambda p, x, y: (fp(p, x) - y)

# using scipy least squares
l1, s =  optimize.leastsq(e, -0.004, args=(x,y))
print l1
# [-0.0132281]


# using numpy least squares
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0]
print l2
# -0.00313461628963 (same answer as Excel trend line)

# smooth x for plotting
x_ = np.arange(0, x[-1], 0.2)

plt.figure()
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-')
plt.show()

编辑 - 附加信息

以上的MWE只包含数据集的一个小样本。当拟合实际数据时，scipy.optimize.curve_fit曲线的R^2值为0.82，而与Excel计算得出的相同的numpy.linalg.lstsq曲线的R^2值为0.41。

- StacyR

Answer 1

您正在最小化不同的误差函数。

当您使用numpy.linalg.lstsq时，被最小化的误差函数为

np.sum((np.log(y) - p * x)**2)

当使用 scipy.optimize.leastsq 函数来最小化一个函数时

np.sum((y - np.exp(p * x))**2)

第一种情况需要因变量和自变量之间的线性依赖关系，但解决方案是已知的分析方法，而第二种情况可以处理任何依赖关系，但依赖于迭代方法。

另外，当使用numpy.linalg.lstsq时，不需要将一行零进行vstack，以下方式同样有效：

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0]

- Jaime

---

谢谢@Jaime - 很棒的回答！不幸的是，我的数学知识不是很好；一个是对的还是错的[也请参见上面的编辑]，或者它们只是根本不同...？对于其他函数有什么影响，例如，如果我想测试将Sigmoid或Gompertz曲线拟合到相同数据的适合性呢？ - StacyR

@StacyR 我没有足够的知识来正确回答你的问题，但我相信像你用 np.linalg.lstsq 拟合指数函数一样只是一个快速而不精确的技巧，无法正确计算误差。这里有一些讨论（对我来说很难理解）：http://mathworld.wolfram.com/LeastSquaresFittingExponential.html 如果你不想深入研究这个东西，我建议使用 scipy 的方法：它应该能够提供更好的拟合效果，并且你的结果将对所有函数保持一致。 - Jaime

再次感谢！我对此进行了更多的研究，正如你所提到的，发现np.linalg.lstsq方法在低x值处过度加权y误差。你分享的链接和我找到的其他资源让我得出了另一种分析方法（使它棘手的是约束条件——所有的书都描述了y=ae^bx的方法而不是y=e^b*x），然而，这也产生了比迭代的scipy.optimize.leastsq更糟糕的拟合曲线。 - StacyR

Answer 2

稍微解释一下Jaime的观点，对数据进行任何非线性转换都会导致不同的误差函数，因此会得到不同的解决方案。这将导致拟合参数的不同置信区间。因此，您有三个可能使用的标准来做出决策：您要最小化哪个误差、您需要更多置信度的参数以及如果您使用拟合来预测某个值，则哪种方法在感兴趣的预测值中产生的误差较小。通过分析和Excel的试验表明，数据中不同类型的噪声（例如，如果噪声函数缩放幅度、影响时间常数或是加性的）会导致选择不同的解决方案。

我还要补充一点，虽然这个技巧对于指数衰减到0的情况“有效”，但它不能用于阻尼指数（上升或下降）通常情况下无法假定为0的值。

- user3117404