首页 文章

SciPy中的指数曲线拟合

提问于
浏览 次
23

我有两个NumPy数组 x 和 y . 当我尝试使用指数函数和 curve_fit （SciPy）使用这个简单的代码来拟合我的数据时

#!/usr/bin/env python
from pylab import *
from scipy.optimize import curve_fit

x = np.array([399.75, 989.25, 1578.75, 2168.25, 2757.75, 3347.25, 3936.75, 4526.25, 5115.75, 5705.25])
y = np.array([109,62,39,13,10,4,2,0,1,2])

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

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

我得错了系数 popt

[a,b,c,d] = [1., 1., 1., 24.19999988]

问题是什么？

python numpy scipy curve-fitting

2 回答

  • 36

    第一条评论：自 a*exp(b - c*x) = (a*exp(b))*exp(-c*x) = A*exp(-c*x) ， a 或 b 是多余的 . 我会放弃 b 并使用：

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

    这不是主要问题 . 问题很简单，当你使用默认的初始猜测（全是1）时， curve_fit 无法收敛到这个问题的解决方案 . 检查 pcov ;你会看到它是 inf . 这并不奇怪，因为如果 c 为1，则 exp(-c*x) 的大部分值下溢为0：

    In [32]: np.exp(-x)
Out[32]: 
array([  2.45912644e-174,   0.00000000e+000,   0.00000000e+000,
         0.00000000e+000,   0.00000000e+000,   0.00000000e+000,
         0.00000000e+000,   0.00000000e+000,   0.00000000e+000,
         0.00000000e+000])

    这表明 c 应该很小 . 更好的初始猜测是，例如 p0 = (1, 1e-6, 1) . 然后我得到：

    In [36]: popt, pcov = curve_fit(func, x, y, p0=(1, 1e-6, 1))

In [37]: popt
Out[37]: array([  1.63561656e+02,   9.71142196e-04,  -1.16854450e+00])

    这看起来很合理：

    In [42]: xx = np.linspace(300, 6000, 1000)

In [43]: yy = func(xx, *popt)

In [44]: plot(x, y, 'ko')
Out[44]: [<matplotlib.lines.Line2D at 0x41c5ad0>]

In [45]: plot(xx, yy)
Out[45]: [<matplotlib.lines.Line2D at 0x41c5c10>]

    fit plot

    回复于
  • 4

    首先，我建议将等式修改为 a*np.exp(-c*(x-b))+d ，否则指数将始终以 x=0 为中心，但并非总是如此 . 您还需要指定合理的初始条件（ curve_fit 的第4个参数指定 [a,b,c,d] 的初始条件） .

    这段代码很合适：

    from pylab import *
from scipy.optimize import curve_fit

x = np.array([399.75, 989.25, 1578.75, 2168.25, 2757.75, 3347.25, 3936.75, 4526.25, 5115.75, 5705.25])
y = np.array([109,62,39,13,10,4,2,0,1,2])

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

popt, pcov = curve_fit(func, x, y, [100,400,0.001,0])
print popt

plot(x,y)
x=linspace(400,6000,10000)
plot(x,func(x,*popt))
show()
    回复于

相关问题