Home Articles

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

Asked
Viewed 1939 times
20

有没有人知道一个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”,它就会失败 .

8 Answers

  • 42

    正确的方法是进行Prony估计,并将结果用作最小二乘拟合的初始猜测(或其他一些更稳健的拟合程序) . Prony估计不需要初始猜测,但确实需要很多点才能产生良好的估计 .

    这是一个概述

    http://www.statsci.org/other/prony.html

    在Octave中,这实现为 expfit ,因此您可以基于Octave库函数编写自己的例程 .

    Prony估计确实需要知道偏移量,但是如果你对衰减“足够远”,你就可以对偏移量进行合理的估计,因此你可以移动数据以将偏移量设置为0.无论如何,Prony估计只是获得其他拟合程序的合理初始猜测的一种方法 .

  • 7

    我会使用 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值 .

  • 6

    我不知道python,但我知道一种简单的方法来非迭代地估计具有偏移的指数衰减系数,给定三个数据点,它们的独立坐标具有固定的差异 . 您的数据点在其独立坐标上有固定的差异(您的x值间隔为60),因此我的方法可以应用于它们 . 你肯定可以将数学翻译成python .

    假设

    y = A + B*exp(-c*x) = A + B*C^x
    

    哪里 C = exp(-c)

    给定y_0,y_1,y_2,对于x = 0,1,2,我们求解

    y_0 = A + B
    y_1 = A + B*C
    y_2 = A + B*C^2
    

    找到A,B,C如下:

    A = (y_0*y_2 - y_1^2)/(y_0 + y_2 - 2*y_1)
    B = (y_1 - y_0)^2/(y_0 + y_2 - 2*y_1)
    C = (y_2 - y_1)/(y_1 - y_0)
    

    相应的指数准确地通过三个点(0,y_0),(1,y_1)和(2,y_2) . 如果您的数据点不在x坐标0,1,2处,而是在k,k s和k 2 * s处,那么

    y = A′ + B′*C′^(k + s*x) = A′ + B′*C′^k*(C′^s)^x = A + B*C^x
    

    所以你可以使用上面的公式来找到A,B,C然后计算

    A′ = A
    C′ = C^(1/s)
    B′ = B/(C′^k)
    

    得到的系数对y坐标中的误差非常敏感,如果推断超出由三个使用的数据点定义的范围,则可能导致大的误差,因此最好从三个数据点计算A,B,C尽可能远(尽管它们之间仍有固定的距离) .

    您的数据集有10个等距数据点 . 让我们选择三个数据点(110,2391),(350,786),(590,263)以供使用 - 它们在独立坐标中具有最大可能的固定距离(240) . 因此,y_0 = 2391,y_1 = 786,y_2 = 263,k = 110,s = 240.然后A = 10.20055,B = 2380.799,C = 0.3258567,A'= 10.20055,B'= 3980.329,C'= 0.9953388 . 指数是

    y = 10.20055 + 3980.329*0.9953388^x = 10.20055 + 3980.329*exp(-0.004672073*x)
    

    您可以将此指数用作非线性拟合算法中的初始猜测 .

    计算A的公式与Shanks变换(http://en.wikipedia.org/wiki/Shanks_transformation)使用的公式相同 .

  • 3

    我从来没有让curve_fit正常工作,因为你说我不想猜测任何东西 . 我试图简化Joe Kington的例子,这就是我的工作 . 我们的想法是将“嘈杂”数据转换为日志,然后将其转换回来并使用polyfit和polyval来确定参数:

    model = np.polyfit(xVals, np.log(yVals) , 1);   
    splineYs = np.exp(np.polyval(model,xVals[0]));
    pyplot.plot(xVals,yVals,','); #show scatter plot of original data
    pyplot.plot(xVals,splineYs('b-'); #show fitted line
    pyplot.show()
    

    其中xVals和yVals只是列表 .

  • 1

    您有两种选择:

    第一种选择是迄今为止最快且最强大的选择 . 但是,它要求您先了解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 才能成为线性系统 .

    如果使用非线性方法,则必须事先知道 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()
    

    Fitting exp

    请注意,线性解决方案提供的结果更接近实际值 . 但是,我们必须提供y偏移值才能使用线性解决方案 . 非线性解决方案不需要这种先验知识 .

  • 1

    @JJacquelin _1841883的答案的Python实现真的很有帮助 . 最初的问题是作为python numpy / scipy请求提出的 . 我把@ johanvdw很干净的R代码并重构为python / numpy . 希望对某人有用:https://gist.github.com/friendtogeoff/00b89fa8d9acc1b2bdf3bdb675178a29

    import numpy as np
    
    """
    compute an exponential decay fit to two vectors of x and y data
    result is in form y = a + b * exp(c*x).
    ref. https://gist.github.com/johanvdw/443a820a7f4ffa7e9f8997481d7ca8b3
    """
    def exp_est(x,y):
        n = np.size(x)
        # sort the data into ascending x order
        y = y[np.argsort(x)]
        x = x[np.argsort(x)]
    
        Sk = np.zeros(n)
    
        for n in range(1,n):
            Sk[n] = Sk[n-1] + (y[n] + y[n-1])*(x[n]-x[n-1])/2
        dx = x - x[0]
        dy = y - y[0]
    
        m1 = np.matrix([[np.sum(dx**2), np.sum(dx*Sk)],
                        [np.sum(dx*Sk), np.sum(Sk**2)]])
        m2 = np.matrix([np.sum(dx*dy), np.sum(dy*Sk)])
    
        [d, c] = (m1.I * m2.T).flat
    
        m3 = np.matrix([[n,                  np.sum(np.exp(  c*x))],
                        [np.sum(np.exp(c*x)),np.sum(np.exp(2*c*x))]])
    
        m4 = np.matrix([np.sum(y), np.sum(y*np.exp(c*x).T)])
    
        [a, b] = (m3.I * m4.T).flat
    
        return [a,b,c]
    
  • 0

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

    enter image description here

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

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

    示例:

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

    enter image description here

    上图是Kington图的副本 .

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

  • 0

    如果你的衰变不是从0开始使用:

    popt, pcov = curve_fit(self.func, x-x0, y)
    

    其中x0是衰变的开始(你想要开始拟合的地方) . 然后再次使用x0进行绘图:

    plt.plot(x, self.func(x-x0, *popt),'--r', label='Fit')
    

    功能是:

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

Related