无初始猜测拟合指数衰减
有没有人知道哪个scipy/numpy模块可以用来把指数衰减拟合到数据上？
我在谷歌上搜索了一下，找到了几个博客文章，比如这个 - http://exnumerus.blogspot.com/2010/04/how-to-fit-exponential-decay-example-in.html，不过这个解决方案需要提前指定y的偏移量，但这并不总是可行的。
编辑：
curve_fit这个函数可以用，但如果没有参数的初始猜测，它可能会失败得很惨，而有时候这个初始猜测是需要的。我正在使用的代码是
#!/usr/bin/env python
import numpy as np
import scipy as sp
import pylab as pl
from scipy.optimize.minpack import curve_fit
x = np.array([ 50., 110., 170., 230., 290., 350., 410., 470.,
530., 590.])
y = np.array([ 3173., 2391., 1726., 1388., 1057., 786., 598.,
443., 339., 263.])
smoothx = np.linspace(x[0], x[-1], 20)
guess_a, guess_b, guess_c = 4000, -0.005, 100
guess = [guess_a, guess_b, guess_c]
exp_decay = lambda x, A, t, y0: A * np.exp(x * t) + y0
params, cov = curve_fit(exp_decay, x, y, p0=guess)
A, t, y0 = params
print "A = %s\nt = %s\ny0 = %s\n" % (A, t, y0)
pl.clf()
best_fit = lambda x: A * np.exp(t * x) + y0
pl.plot(x, y, 'b.')
pl.plot(smoothx, best_fit(smoothx), 'r-')
pl.show()
这个代码是可以工作的，但如果我们去掉“p0=guess”，它就会失败得很惨。
9 个回答
我会使用 scipy.optimize.curve_fit 这个函数。它的说明文档里甚至有一个关于如何拟合指数衰减的例子，我在这里复制过来：
>>> import numpy as np
>>> 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)
因为加了随机噪声，所以拟合出来的参数会有所不同，但我得到了 2.47990495、1.40709306 和 0.53753635 作为 a、b 和 c，这在有噪声的情况下还算不错。如果我用 y 来拟合，而不是 yn，我就能得到准确的 a、b 和 c 值。
没有初始猜测的指数拟合过程，不是迭代过程：
这段内容来自于论文（第16-17页）：https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales
如果需要，这可以用来初始化一个非线性回归计算，以选择特定的优化标准。
示例：
Joe Kington给出的例子很有趣。不幸的是，数据没有显示，只有图表。因此，下面的（x,y）数据来自于对图表的扫描，所以这些数值可能和Joe Kington使用的数值不完全一样。不过，"拟合"曲线的方程非常接近，考虑到数据点的分散程度。
上面的图是Kington的图表的复制。
下面的图展示了使用上述过程得到的结果。
更新：一种变体
你有两个选择：
- 把系统线性化，然后对数据的对数进行线性拟合。
- 使用非线性求解器（比如
scipy.optimize.curve_fit）。
第一个选择是最快且最可靠的。不过，它需要你事先知道y的偏移量，否则就无法把方程线性化。（也就是说，y = A * exp(K * t)可以通过拟合y = log(A * exp(K * t)) = K * t + log(A)来线性化，但y = A*exp(K*t) + C只能通过拟合y - C = K*t + log(A)来线性化，而因为y是你的自变量，所以C必须事先知道，这样才能形成线性系统。）
如果你使用非线性的方法，a) 不保证能收敛并找到解，b) 会慢很多，c) 对参数的不确定性估计会差很多，d) 通常精度也会低很多。不过，非线性方法有一个很大的优点：它可以解决非线性方程组。在你的情况下，这意味着你不需要事先知道C。
举个例子，我们用一些带噪声的数据来求解y = A * exp(K * t)，分别使用线性和非线性的方法：
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.optimize
def main():
# Actual parameters
A0, K0, C0 = 2.5, -4.0, 2.0
# Generate some data based on these
tmin, tmax = 0, 0.5
num = 20
t = np.linspace(tmin, tmax, num)
y = model_func(t, A0, K0, C0)
# Add noise
noisy_y = y + 0.5 * (np.random.random(num) - 0.5)
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
# Non-linear Fit
A, K, C = fit_exp_nonlinear(t, noisy_y)
fit_y = model_func(t, A, K, C)
plot(ax1, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, C0))
ax1.set_title('Non-linear Fit')
# Linear Fit (Note that we have to provide the y-offset ("C") value!!
A, K = fit_exp_linear(t, y, C0)
fit_y = model_func(t, A, K, C0)
plot(ax2, t, y, noisy_y, fit_y, (A0, K0, C0), (A, K, 0))
ax2.set_title('Linear Fit')
plt.show()
def model_func(t, A, K, C):
return A * np.exp(K * t) + C
def fit_exp_linear(t, y, C=0):
y = y - C
y = np.log(y)
K, A_log = np.polyfit(t, y, 1)
A = np.exp(A_log)
return A, K
def fit_exp_nonlinear(t, y):
opt_parms, parm_cov = sp.optimize.curve_fit(model_func, t, y, maxfev=1000)
A, K, C = opt_parms
return A, K, C
def plot(ax, t, y, noisy_y, fit_y, orig_parms, fit_parms):
A0, K0, C0 = orig_parms
A, K, C = fit_parms
ax.plot(t, y, 'k--',
label='Actual Function:\n $y = %0.2f e^{%0.2f t} + %0.2f$' % (A0, K0, C0))
ax.plot(t, fit_y, 'b-',
label='Fitted Function:\n $y = %0.2f e^{%0.2f t} + %0.2f$' % (A, K, C))
ax.plot(t, noisy_y, 'ro')
ax.legend(bbox_to_anchor=(1.05, 1.1), fancybox=True, shadow=True)
if __name__ == '__main__':
main()
注意，线性解的结果与实际值更接近。不过，我们必须提供y的偏移量值才能使用线性解。非线性解则不需要这些事先的知识。