在SciPy中进行指数曲线拟合

Question

我有两个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]

什么是问题？

- drastega

---

类似的问题：http://stackoverflow.com/questions/17527869/curve-fit-fails-with-exponential-but-zunzun-gets-it-right - Josef

Answer 1

首先评论：由于a*exp(b - c*x) = (a*exp(b))*exp(-c*x) = A*exp(-c*x)，所以a或b是多余的。我会去掉b并使用：

import matplotlib.pyplot as plt

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]: plt.plot(x, y, 'ko')
Out[44]: [<matplotlib.lines.Line2D at 0x41c5ad0>]

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

- Warren Weckesser

---

为什么要使用-c而不是c？如果需要，curve_fit可以找到负c，不是吗？ - Elliot Gorokhovsky

@RenéG：这是drastega在问题中使用的约定。 - Warren Weckesser

另一种处理初始参数的方法（使用默认值）是将_x_归一化为（大约）0-1，例如，_ξ=x/k_，估计_a_、_c'_和_d_，最终得到_c=c'/k_。 - gboffi

Answer 2

首先我建议你修改你的方程为a*np.exp(-c*(x-b))+d，否则指数函数总是以x=0为中心，这可能并不总是正确的。

还需要指定合理的初始条件（curve_fit的第四个参数指定了[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()

- three_pineapples

---

初始条件从何而来？ - Marcin Zdunek

@MarcinZdunek 这是一段时间以前的事情，所以我记不清楚了。振幅将会从图形中被估计出来。其它的可能会通过试错法被确定，尽管 c 的值也可以被估算出来（请参考这个问题的被接受答案）。 - three_pineapples

@MarcinZdunek 如果您对两个数据范围进行归一化，然后再对估计参数进行反归一化，则默认初始值是可以的... - gboffi

我想强调一下，重新查看一遍后，我认为 a 和 b 的初始条件来自第一个 y 和 x 值（假设值是按顺序排列的），可以像接受的答案中那样估计 c ，而对于 d 的估计来自最终的 y 值，这些值约为 ~0。如果你对初始条件感到困惑，这可以成为一个良好的起点。 - three_pineapples