在没有初始猜测的情况下拟合指数衰减

Question

有没有人知道一个scipy/numpy模块,它将允许指数衰减数据？

谷歌搜索返回了一些博客文章,例如http://exnumerus.blogspot.com/2010/04/how-to-fit-exponential-decay-example-in.html,但该解决方案要求y-offset为预先指定的,这并不总是可行的

编辑:

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",它就会失败.

Answer 1

您有两种选择:

1. 线性化系统,并在数据日志中插入一条线.
2. 使用非线性求解器(例如 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偏移值才能使用线性解决方案.非线性解决方案不需要这种先验知识.

Answer 2

我会用这个scipy.optimize.curve_fit功能.它的doc字符串甚至有一个在其中拟合指数衰减的例子,我将在这里复制:

>>> 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值.

Answer 3

在没有初始猜测而不是迭代过程的情况下拟合指数的过程:

这来自论文(第16-17页):https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

如有必要,可以使用此方法初始化非线性回归演算,以便选择特定的优化标准.

示例:

Joe Kington给出的例子很有趣.不幸的是,数据没有显示,只有图表.因此,下面的数据(x,y)来自图形的图形扫描,因此数值可能不是Joe Kington使用的数值.然而,考虑到点的广泛分散,"拟合"曲线的各个方程彼此非常接近.

上图是Kington图的副本.

下图显示了使用上述程序获得的结果.