在没有初始猜测的情况下拟合指数衰减-Java 学习之路

有没有人知道一个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 回答

您有两种选择：

线性化系统，并在数据日志中插入一行 .
使用非线性求解器（例如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 才能成为线性系统 .

如果使用非线性方法，则必须事先知道 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偏移值才能使用线性解决方案 . 非线性解决方案不需要这种先验知识 .

回复于 2024-04-20T23:32:58+08:00

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值 .
回复于 2024-04-20T23:32:58+08:00
0

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

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

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

示例：

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

上图是Kington图的副本 .

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

回复于 2024-04-20T23:32:58+08:00
42

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

这是一个概述

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

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

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

回复于 2024-04-20T23:32:58+08:00
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只是列表 .
回复于 2024-04-20T23:32:58+08:00
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）使用的公式相同 .
回复于 2024-04-20T23:32:58+08:00
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
```
回复于 2024-04-20T23:32:58+08:00

@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]

回复于 2024-04-20T23:32:58+08:00

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

8 回答

相关问题