无需初始猜测即可拟合指数衰减 [英] fitting exponential decay with no initial guessing
问题描述
有谁知道一个 scipy/numpy 模块可以让数据适应指数衰减?
Does anyone know a scipy/numpy module which will allow to fit exponential decay to data?
Google 搜索返回了一些博客文章，例如 - http://exnumerus.blogspot.com/2010/04/how-to-fit-exponential-decay-example-in.html ，但该解决方案需要预先设置 y-offset指定，这并不总是可能的
Google search returned a few blog posts, for example - http://exnumerus.blogspot.com/2010/04/how-to-fit-exponential-decay-example-in.html , but that solution requires y-offset to be pre-specified, which is not always possible
curve_fit 有效，但它可能会在没有对参数的初始猜测的情况下非常悲惨地失败，而这有时是需要的.我正在使用的代码是
curve_fit works, but it can fail quite miserably with no initial guess for parameters, and that is sometimes needed. The code I'm working with is
#!/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
t = %s
y0 = %s
" % (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"，它就会惨遭失败.
which works, but if we remove "p0=guess", it fails miserably.
推荐答案
您有两个选择:
- 使系统线性化，并在数据日志中拟合一条线.
- 使用非线性求解器(例如
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 才能使其成为线性系统.
The first option is by far the fastest and most robust. However, it requires that you know the y-offset a-priori, otherwise it's impossible to linearize the equation. (i.e. y = A * exp(K * t) can be linearized by fitting y = log(A * exp(K * t)) = K * t + log(A), but y = A*exp(K*t) + C can only be linearized by fitting y - C = K*t + log(A), and as y is your independent variable, C must be known beforehand for this to be a linear system.
如果您使用非线性方法，则 a) 不能保证收敛并产生解决方案，b) 会慢得多，c) 对参数中的不确定性的估计要差得多，d) 通常是不那么精确.然而，非线性方法与线性反演相比有一个巨大的优势:它可以求解非线性方程组.就您而言，这意味着您不必事先了解 C.
If you use a non-linear method, it's a) not guaranteed to converge and yield a solution, b) will be much slower, c) gives a much poorer estimate of the uncertainty in your parameters, and d) is often much less precise. However, a non-linear method has one huge advantage over a linear inversion: It can solve a non-linear system of equations. In your case, this means that you don't have to know C beforehand.
举个例子，让我们使用线性和非线性方法用一些噪声数据求解 y = A * exp(K * t):
Just to give an example, let's solve for y = A * exp(K * t) with some noisy data using both linear and nonlinear methods:
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:
$y = %0.2f e^{%0.2f t} + %0.2f$' % (A0, K0, C0))
ax.plot(t, fit_y, 'b-',
label='Fitted Function:
$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 偏移值才能使用线性解决方案.非线性解决方案不需要这种先验知识.
Note that the linear solution provides a result much closer to the actual values. However, we have to provide the y-offset value in order to use a linear solution. The non-linear solution doesn't require this a-priori knowledge.
这篇关于无需初始猜测即可拟合指数衰减的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
