python - 如何在 Python 中进行指数和对数曲线拟合？我发现只有多项式拟合

Question

我有一组数据，我想比较哪一行最能描述它(不同阶的多项式、指数或对数)。

我使用 Python 和 Numpy，对于多项式拟合，有一个函数 polyfit()。但是我发现没有这样的函数用于指数和对数拟合。

有吗？或者如何解决？

Answer 1

为了拟合y = A + B log x，正好拟合y 反对 (log x)。

>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> numpy.polyfit(numpy.log(x), y, 1)
array([ 8.46295607,  6.61867463])
# y ≈ 8.46 log(x) + 6.62

---

为了拟合y = Ae^Bx，取两边的对数给出log y = log A + Bx。所以适合 (log y) 对 x。

请注意，将 (log y) 拟合为线性会强调 y 的较小值，从而导致较大的 y 偏差较大。这是因为 polyfit(线性回归)通过最小化 ∑_i (ΔY)² = ∑_i (Y_i - Ŷ_{< em>i})²。当Y_i = log y_i，残基 ΔY_i = Δ(log y_i) ≈ Δy_i/|y_i |。因此，即使 polyfit 对较大的 y 做出了非常糟糕的决定，“除以 -|y|” factor 将对其进行补偿，导致 polyfit 偏爱小值。

这可以通过给每个条目一个与y成比例的“权重”来缓解。 polyfit 通过 w 关键字参数支持加权最小二乘。

>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> numpy.polyfit(x, numpy.log(y), 1)
array([ 0.10502711, -0.40116352])
#    y ≈ exp(-0.401) * exp(0.105 * x) = 0.670 * exp(0.105 * x)
# (^ biased towards small values)
>>> numpy.polyfit(x, numpy.log(y), 1, w=numpy.sqrt(y))
array([ 0.06009446,  1.41648096])
#    y ≈ exp(1.42) * exp(0.0601 * x) = 4.12 * exp(0.0601 * x)
# (^ not so biased)

请注意，Excel、LibreOffice 和大多数科学计算器通常对指数回归/趋势线使用未加权(有偏)公式。如果您希望您的结果与这些平台兼容，请不要包括权重，即使它提供了更好的结果。

---

现在，如果你可以使用 scipy，你可以使用 scipy.optimize.curve_fit无需转换即可拟合任何模型。

对于y = A + B log x 结果与变换方法相同:

>>> x = numpy.array([1, 7, 20, 50, 79])
>>> y = numpy.array([10, 19, 30, 35, 51])
>>> scipy.optimize.curve_fit(lambda t,a,b: a+b*numpy.log(t),  x,  y)
(array([ 6.61867467,  8.46295606]), 
 array([[ 28.15948002,  -7.89609542],
        [ -7.89609542,   2.9857172 ]]))
# y ≈ 6.62 + 8.46 log(x)

对于 y = Ae^Bx，我们可以得到更好的拟合，因为它计算 Δ(直接记录y)。但是我们需要提供一个初始化猜测，以便 curve_fit 可以达到所需的局部最小值。

>>> x = numpy.array([10, 19, 30, 35, 51])
>>> y = numpy.array([1, 7, 20, 50, 79])
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y)
(array([  5.60728326e-21,   9.99993501e-01]),
 array([[  4.14809412e-27,  -1.45078961e-08],
        [ -1.45078961e-08,   5.07411462e+10]]))
# oops, definitely wrong.
>>> scipy.optimize.curve_fit(lambda t,a,b: a*numpy.exp(b*t),  x,  y,  p0=(4, 0.1))
(array([ 4.88003249,  0.05531256]),
 array([[  1.01261314e+01,  -4.31940132e-02],
        [ -4.31940132e-02,   1.91188656e-04]]))
# y ≈ 4.88 exp(0.0553 x). much better.