如何使用python拟合多条指数曲线

问题描述

对于单个指数曲线(如此处的图片curve_fit for as single exponential curve所示),我可以使用scipy.optimize.curve_fit拟合数据。但是,我不确定如何实现由多条指数曲线组成的相似数据集的拟合,如此处double exponential curves所示。 我使用以下方法实现了单曲线的拟合:

def exp_decay(x,a,r):
    return a * ((1-r)**x) 

x = np.linspace(0,50,50)
y = exp_decay(x,400,0.06)

y1 = exp_decay(x,550,0.06)      # this is to be used to append to y to generate two curves

pars,cov = curve_fit(exp_decay,x,y,p0=[0,0])
plt.scatter(x,y)
plt.plot(x,exp_decay(x,*pars),'r-')     #this realizes the fit for a single curve

yx = np.append(y,y1)   #this realizes two exponential curves (as shown above - double exponential curves) for which I don't need to fit a model to

有人可以帮助描述如何针对两条曲线的数据集实现这一目标。我的实际数据集包含多条指数曲线,但我认为,如果我可以实现两条曲线的拟合,则可以为我的数据集复制相同的曲线。绝不能使用scipy的curve_fit来完成;任何可行的实现方式都可以。

请帮助!!!

解决方法

您可以通过使用简单的准则(例如一阶导数估计)分割数据集来轻松解决您的问题,然后我们可以对每个子数据集应用简单的曲线拟合过程。

试验数据集

首先,让我们导入一些程序包并创建具有三个曲线的综合数据集来表示您的问题。

我们使用两个参数的指数模型,因为时间原点偏移将由拆分方法处理。我们还会添加噪音,因为现实世界的数据总是有噪音:

import numpy as np
import pandas as pd
from scipy import optimize
import matplotlib.pyplot as plt

def func(x,a,b):
    return a*np.exp(b*x)

N = 1001
n1 = N//3
n2 = 2*n1

t = np.linspace(0,10,N)

x0 = func(t[:n1],1,-0.2)
x1 = func(t[n1:n2]-t[n1],5,-0.4)
x2 = func(t[n2:]-t[n2],2,-1.2)

x = np.hstack([x0,x1,x2])
xr = x + 0.025*np.random.randn(x.size)

以图形方式呈现如下:

enter image description here

数据集拆分

我们可以使用一个简单的准则将数据集分为三个子数据集,作为使用一阶差异进行评估的一阶导数估计。目的是检测曲线何时急剧上升或下降(应该分割数据集。估计一阶导数如下):

dxrdt = np.abs(np.diff(xr)/np.diff(t))

该标准需要一个额外的参数(阈值),必须根据您的信号规格进行相应的调整。该标准等效于:

xcrit = 20
q = np.where(dxrdt > xcrit) # (array([332,665],dtype=int64),)

和分割索引为:

idx = [0] + list(q[0]+1) + [t.size] # [0,333,666,1001]

主要是,标准阈值将受到数据中噪声的性质和功率以及两条曲线之间的间隙幅度的影响。该方法的使用取决于在存在噪声的情况下检测曲线间隙的能力。当噪声功率与我们要检测的间隙大小相同时,它将破裂。如果噪声严重拖尾(极少的异常值),您也可以观察到虚假的分割指数。

在此MCVE中,我们将阈值设置为20 [Signal Units/Time Units]

enter image description here

对此手工制定的标准的另一种选择是将标识委托给scipy的出色find_peaks方法。但这并不能避免将检测调整到信号规格的要求。

适合原点偏移的数据集

现在,我们可以将曲线拟合应用于每个子数据集(原点偏移时间),收集参数和统计数据并绘制结果:

trials = []
fig,axe = plt.subplots()
for k,(i,j) in enumerate(zip(idx[:-1],idx[1:])):
    p,s = optimize.curve_fit(func,t[i:j]-t[i],xr[i:j])
    axe.plot(t[i:j],xr[i:j],'.',label="Data #{}".format(k+1))
    axe.plot(t[i:j],func(t[i:j]-t[i],*p),label="Data Fit #{}".format(k+1))
    trials.append({"n0": i,"n1": j,"t0": t[i],"a": p[0],"b": p[1],"s_a": s[0,0],"s_b": s[1,1],"s_ab": s[0,1]})
axe.set_title("Curve Fits")
axe.set_xlabel("Time,$t$")
axe.set_ylabel("Signal Estimate,$\hat{g}(t)$")
axe.legend()
axe.grid()
df = pd.DataFrame(trials)

它返回以下拟合结果:

    n0    n1    t0         a         b       s_a           s_b      s_ab
0    0   333  0.00  0.998032 -0.199102  0.000011  4.199937e-06 -0.000005
1  333   666  3.33  5.001710 -0.399537  0.000013  3.072542e-07 -0.000002
2  666  1001  6.66  2.002495 -1.203943  0.000030  2.256274e-05 -0.000018

符合我们的原始参数(请参见“试验数据集”部分)。

以图形方式我们可以检查拟合的优度:

enter image description here