python函数拟合不规则曲线_如何在Python中做指数和对数曲线拟合?我发现只有多项式拟合...

我有一组数据,我想比较哪一种曲线可以最好地描述它(不同阶的多项式,指数或对数)。

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

有没有?或者如何解决呢?

最佳解决思路

为了拟合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 = AeBx,取双方的对数给出logy = logA + Bx。所以适合(logy)对x。

请注意,拟合(log y)好像是线性的,会强调y的小值,导致大y的偏差较大。这是因为polyfit(线性回归)通过最小化Σi(ΔY)2 =Σi(Yi-Ŷi)2来工作。当Yi = log yi时,残差ΔYi=Δ(log yi)≈Δyi/| yi |。所以即使polyfit对于大y做出了一个非常糟糕的决定,"divide-by-|y|"因子也会补偿它,导致polyfit偏好较小的值。

这可以通过给每个条目一个与y成比例的"weight"来减轻。 polyfit通过w关键字参数支持weighted-least-squares。

>>> 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 = AeBx,我们可以更好地拟合,因为它直接计算Δ(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.

b08b2ad9e230bdeb81d214a819fc0829.jpg

次佳解决思路

您也可以使用来自scipy.optimize的curve_fit来拟合一组数据。例如,如果你想拟合一个指数函数(来自documentation):

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)

然后如果你想要绘制图像(plot),你可以这样做:

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()

(注意:当你绘图的时候,*之前的*会将这些条件扩展到func所期望的a,b和c中。)

第三种解决思路

我在这方面遇到了一些麻烦,所以让我说得更细致清楚一些,以便像我这样的新手可以理解。

比方说,我们有一个数据文件或类似的东西

# -*- 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

beaef97daf51559cc4cae6f2b29306da.jpg

参考资料

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值