如何在Python中进行指数和对数曲线拟合？我只找到了多项式拟合

Question

我有一组数据，想要比较哪一条线最能准确描述这些数据（比如不同阶数的多项式、指数函数或者对数函数）。

我在用Python和Numpy，针对多项式拟合有一个叫polyfit()的函数。但是我没有找到适合指数和对数拟合的函数。

有没有这样的函数呢？或者还有其他解决办法吗？

Answer 1

我在这个问题上遇到了一些麻烦，所以我想把事情说得很清楚，让像我这样的菜鸟也能理解。

假设我们有一个数据文件或者类似的东西。

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import sympy as sym

"""
Generate some data, let's imagine that you already have this. 
"""
x = np.linspace(0, 3, 50)
y = np.exp(x)

"""
Plot your data
"""
plt.plot(x, y, 'ro',label="Original Data")

"""
brutal force to avoid errors
"""    
x = np.array(x, dtype=float) #transform your data in a numpy array of floats 
y = np.array(y, dtype=float) #so the curve_fit can work

"""
create a function to fit with your data. a, b, c and d are the coefficients
that curve_fit will calculate for you. 
In this part you need to guess and/or use mathematical knowledge to find
a function that resembles your data
"""
def func(x, a, b, c, d):
    return a*x**3 + b*x**2 +c*x + d

"""
make the curve_fit
"""
popt, pcov = curve_fit(func, x, y)

"""
The result is:
popt[0] = a , popt[1] = b, popt[2] = c and popt[3] = d of the function,
so f(x) = popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3].
"""
print "a = %s , b = %s, c = %s, d = %s" % (popt[0], popt[1], popt[2], popt[3])

"""
Use sympy to generate the latex sintax of the function
"""
xs = sym.Symbol('\lambda')    
tex = sym.latex(func(xs,*popt)).replace('$', '')
plt.title(r'$f(\lambda)= %s$' %(tex),fontsize=16)

"""
Print the coefficients and plot the funcion.
"""

plt.plot(x, func(x, *popt), label="Fitted Curve") #same as line above \/
#plt.plot(x, popt[0]*x**3 + popt[1]*x**2 + popt[2]*x + popt[3], label="Fitted Curve") 

plt.legend(loc='upper left')
plt.show()

结果是：a = 0.849195983017 , b = -1.18101681765, c = 2.24061176543, d = 0.816643894816

Answer 2

你可以使用 scipy.optimize 里的 curve_fit 来把一组数据拟合成你喜欢的任何函数。比如说，如果你想拟合一个指数函数（具体可以参考文档）：

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

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

x = np.linspace(0,4,50)
y = func(x, 2.5, 1.3, 0.5)
yn = y + 0.2*np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn)

然后如果你想画图，可以这样做：

plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()

（注意：在你画图时，popt 前面的 * 是用来把里面的内容展开成 a、b 和 c，这些是 func 需要的参数。）

Answer 3

要拟合 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 − Ŷ_i)² 来工作的。当 Y_i = log y_i 时，残差 ΔY_i = Δ(log y_i) ≈ Δy_i / |y_i|。所以即使 polyfit 对大的 y 值做了很糟糕的决策，"除以 |y|" 的因素也会弥补这个问题，导致 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，我们可以得到更好的拟合，因为它直接计算 Δ(log 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.