对于单个指数曲线,如下图所示 单个指数曲线的curve_fit,我可以使用scipy.optimize.curve_fit来拟合数据。然而,我不确定如何对类似数据集进行拟合,该数据集由多个指数曲线组成,如下图所示 双指数曲线。我使用以下方法实现了单曲线的拟合:
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) # 这将用于附加到y以生成两条曲线pars, cov = curve_fit(exp_decay, x, y, p0=[0,0])plt.scatter(x,y)plt.plot(x, exp_decay(x, *pars), 'r-') # 这实现了单曲线的拟合yx = np.append(y,y1) # 这实现了两个指数曲线(如上所示 - 双指数曲线),我无需为其拟合模型
能否有人帮助描述如何对包含两条曲线的数据集实现这一点。我的实际数据集包含多个指数曲线,但我认为如果我能实现对两条曲线的拟合,我可能能够对我的数据集进行同样的操作。这不能使用scipy的curve_fit来完成;任何有效的实现都可以。
请帮助我!!!
回答:
你的问题可以通过使用简单的判据(如一阶导数估计)来分割你的数据集,然后对每个子数据集应用简单的曲线拟合程序来轻松解决。
试验数据集
首先,让我们导入一些包并创建一个包含三条曲线的合成数据集来表示你的问题。
我们使用一个两参数指数模型,因为时间原点偏移将由分割方法处理。我们还添加了噪声,因为现实世界的数据总是有噪声的:
import numpy as npimport pandas as pdfrom scipy import optimizeimport matplotlib.pyplot as pltdef func(x, a, b): return a*np.exp(b*x)N = 1001n1 = N//3n2 = 2*n1t = 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)
图形上呈现如下:
数据集分割
我们可以使用一阶导数估计作为简单的判据,通过一阶差分来评估它,将数据集分割成三个子数据集。目标是检测曲线何时急剧上升或下降(即数据集应分割的位置)。一阶导数估计如下:
dxrdt = np.abs(np.diff(xr)/np.diff(t))
判据需要一个额外的参数(阈值),必须根据你的信号规格进行调整。判据相当于:
xcrit = 20q = np.where(dxrdt > xcrit) # (array([332, 665], dtype=int64),)
分割索引为:
idx = [0] + list(q[0]+1) + [t.size] # [0, 333, 666, 1001]
主要是判据阈值将受到数据噪声的性质和强度以及两条曲线之间间隙大小的影响。使用这种方法取决于在存在噪声的情况下检测曲线间隙的能力。当噪声强度与我们想要检测的间隙大小相同的时候,它将失效。如果噪声具有重尾分布(少数强烈的异常值),你也可能会观察到错误的分割索引。
在这个最小完整示例中,我们将阈值设置为 20 [信号单位/时间单位]
:
这种手工制作的判据的替代方案是将识别工作委托给 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="数据 #{}".format(k+1)) axe.plot(t[i:j], func(t[i:j]-t[i], *p), label="数据拟合 #{}".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("曲线拟合")axe.set_xlabel("时间, $t$")axe.set_ylabel("信号估计, $\hat{g}(t)$")axe.legend()axe.grid()df = pd.DataFrame(trials)
它返回以下拟合结果:
n0 n1 t0 a b s_a s_b s_ab0 0 333 0.00 0.998032 -0.199102 0.000011 4.199937e-06 -0.0000051 333 666 3.33 5.001710 -0.399537 0.000013 3.072542e-07 -0.0000022 666 1001 6.66 2.002495 -1.203943 0.000030 2.256274e-05 -0.000018
这与我们原始的参数相符(见试验数据集部分)。
图形上我们可以检查拟合的优劣: