python - 使用 scipy curve_fit 拟合噪声指数的建议？

Question

我正在尝试拟合通常按以下方式建模的数据:

def fit_eq(x, a, b, c, d, e):
    return a*(1-np.exp(-x/b))*(c*np.exp(-x/d)) + e

x = np.arange(0, 100, 0.001)
y = fit_eq(x, 1, 1, -1, 10, 0)
plt.plot(x, y, 'b')

不过，实际轨迹的一个例子要嘈杂得多:

如果我分别拟合上升分量和衰减分量，我可以得出比较合适的结果:

def fit_decay(df, peak_ix):
    fit_sub = df.loc[peak_ix:]

    guess = np.array([-1, 1e-3, 0])
    x_zeroed = fit_sub.time - fit_sub.time.values[0]

    def exp_decay(x, a, b, c):
        return a*np.exp(-x/b) + c

    popt, pcov = curve_fit(exp_decay, x_zeroed, fit_sub.primary, guess)

    fit = exp_decay(x_full_zeroed, *popt)

    return x_zeroed, fit_sub.primary, fit

def fit_rise(df, peak_ix):
        fit_sub = df.loc[:peak_ix]
        guess = np.array([1, 1, 0])
        def exp_rise(x, a, b, c):
             return a*(1-np.exp(-x/b)) + c

        popt, pcov = curve_fit(exp_rise, fit_sub.time, 
                              fit_sub.primary, guess, maxfev=1000)

        x = df.time[:peak_ix+1]
        y = df.primary[:peak_ix+1]
        fit = exp_rise(x.values, *popt)

        return x, y, fit

ix = df.primary.idxmin()

rise_x, rise_y, rise_fit = fit_rise(df, ix)
decay_x, decay_y, decay_fit = fit_decay(df, ix)
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.plot(rise_x, rise_y)
ax1.plot(rise_x, rise_fit)
ax2.plot(decay_x, decay_y)
ax2.plot(decay_x, decay_fit)

不过，理想情况下，我应该能够使用上面的等式拟合整个 transient 。不幸的是，这不起作用:

def fit_eq(x, a, b, c, d, e):
    return a*(1-np.exp(-x/b))*(c*np.exp(-x/d)) + e

guess = [1, 1, -1, 1, 0]

x = df.time
y = df.primary

popt, pcov = curve_fit(fit_eq, x, y, guess)
fit = fit_eq(x, *popt)
plt.plot(x, y)
plt.plot(x, fit)

我已经为 guess 尝试了许多不同的组合，包括我认为应该是合理近似值的数字，但要么我得到了糟糕的拟合，要么 curve_fit 找不到参数。

我也尝试过拟合较小的数据部分(例如 0.12 到 0.16 秒)，但没有取得更大的成功。

此特定示例的数据集副本位于此处 Share CSV

我在这里缺少任何提示或技巧吗？

编辑 1:

因此，如建议的那样，如果我限制适合的区域不包括左侧的高原(即下图中的橙色)，我会得到一个不错的适合。我看到了另一篇关于 curve_fit 的 stackoverflow 帖子，其中提到转换非常小的值也有帮助。将时间变量从秒转换为毫秒对获得合适的拟合产生了很大的影响。

我还发现，强制 curve_fit 尝试通过一些点(特别是峰值，然后是衰减拐点处的一些较大点，因为那里的各种瞬变会拉低衰减拟合) .

我想对于左边的高原，我可以拟合一条线并将其连接到指数拟合？我最终想要实现的是减去大的 transient ，所以我需要在左侧表示一些高原。

sub = df[(df.time>0.1275) & (d.timfe < 0.6)]

def fit_eq(x, a, b, c, d, e):
    return a*(1-np.exp(-x/b))*(np.exp(-x/c) + np.exp(-x/d)) + e 

x = sub.time
x = sub.time - sub.time.iloc[0]
x *= 1e3
y = sub.primary

guess = [-1, 1, 1, 1, -60]
ixs = y.reset_index(drop=True)[100:300].sort_values(ascending=False).index.values[:10]
ixmin = y.reset_index(drop=True).idxmin()
sigma = np.ones(len(x))
sigma[ixs] = 0.1
sigma[ixmin] = 0.1

popt, pcov = curve_fit(fit_eq, x, y, p0=guess, sigma=sigma, maxfev=2000)

fit = fit_eq(x, *popt)
x = x*1e-3

f, (ax1, ax2) = plt.subplots(1,2, figsize=(16,8))
ax1.plot((df.time-sub.time.iloc[0]), df.primary)
ax1.plot(x, y)
ax1.plot(x.iloc[ixs], y.iloc[ixs], 'o')
ax1.plot(x, fit, lw=4)
ax2.plot((df.time-sub.time.iloc[0]), df.primary)
ax2.plot(x, y)
ax2.plot(x.iloc[ixs], y.iloc[ixs], 'o')
ax2.plot(x, fit)
ax1.set_xlim(-.02, .06)

Answer 1

我来自 EE 背景，寻找“系统识别”工具，但没有在我找到的 Python 库中找到我期望的东西

所以我在我比较熟悉的频域中制定了一个“朴素”的SysID解决方案

我删除了初始偏移量，假设阶跃激励，加倍，将数据集反转为 fft 处理步骤的周期性

使用 scipy.optimize.least_squares 拟合拉普拉斯/频域传递函数后:

def tf_model(w, td0,ta,tb,tc): # frequency domain transfer function w delay
    return np.exp(-1j*w/td0)*(1j*w*ta)/(1j*w*tb + 1)/(1j*w*tc + 1)

在 sympy 的帮助下，我转换回了时域阶跃响应

inverse_laplace_transform(s*a/((s*b + 1)*(s*c + 1)*s), s, t

稍微简化后:

def tdm(t, a, b, c):
    return -a*(np.exp(-t/c) - np.exp(-t/b))/(b - c)

对频域拟合常数应用归一化，排列图

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import least_squares

data = np.loadtxt(open("D:\Downloads\\transient_data.csv","rb"),
                  delimiter=",", skiprows=1)

x, y = zip(*data[1:]) # unpacking, dropping one point to get 1000 
x, y = np.array(x), np.array(y)

y = y - np.mean(y[:20]) # remove linear baseline from starting data estimate

xstep = np.sign((x - .12))*-50 # eyeball estimate step start time, amplitude

x = np.concatenate((x,x + x[-1]-x[0])) # extend, invert for a periodic data set
y = np.concatenate((y, -y))
xstep = np.concatenate((xstep, -xstep))

# frequency domain transforms of the data, assumed square wave stimulus
fy = np.fft.rfft(y)
fsq = np.fft.rfft(xstep)

# only keep 1st ~50 components of the square wave
# this is equivalent to applying a rectangular window low pass
K = np.arange(1,100,2) # 1st 50 nonzero fft frequency bins of the square wave
# form the frequency domain transfer function from fft data: Gd
Gd = fy[1:100:2]/fsq[1:100:2]

def tf_model(w, td0,ta,tb,tc): # frequency domain transfer function w delay
    return np.exp(-1j*w/td0)*(1j*w*ta)/(1j*w*tb + 1)/(1j*w*tc + 1)

td0,ta,tb,tc = 0.1, -1, 0.1, 0.01

x_guess = [td0,ta,tb,tc]

# cost function, "residual" with weighting by stimulus frequency components**2?
def func(x, Gd, K):
    return (np.conj(Gd - tf_model(K, *x))*
                   (Gd - tf_model(K, *x))).real/K #/K # weighting by K powers

res = least_squares(func, x_guess, args=(Gd, K),
                    bounds=([0.0, -100, 0, 0],
                            [1.0, 0.0, 10, 1]),
                             max_nfev=100000, verbose=1)

td0,ta,tb,tc = res['x']

# convolve model w square wave in frequency domain
fy = fsq * tf_model(np.arange(len(fsq)), td0,ta,tb,tc)

ym = np.fft.irfft(fy) # back to time domain 

print(res)

plt.plot(x, xstep, 'r')
plt.plot(x, y, 'g')
plt.plot(x, ym, 'k')

# finally show time domain step response function, normaliztion
def tdm(t, a, b, c):
    return -a*(np.exp(-t/c) - np.exp(-t/b))/(b - c)

# normalizing factor for frequency domain, dataset time range
tn = 2*np.pi/(x[-1]-x[0])
ta, tb, tc = ta/tn, tb/tn, tc/tn

y_tdm = tdm(x - 0.1, ta, tb, tc)

# roll shifts yellow y_tdm to (almost) match black frequency domain model
plt.plot(x, 100*np.roll(y_tdm, 250), 'y')

绿色:翻倍，倒置的数据是周期性的
红色:估计的起始步骤，也加倍，反转为周期性方波
黑色:与方波卷积的频域拟合模型
黄色:拟合频域模型转换回时域阶跃响应，滚动比较

     message: '`ftol` termination condition is satisfied.'
        nfev: 40
        njev: 36
  optimality: 0.0001517727368912258
      status: 2
     success: True
           x: array([ 0.10390021, -0.4761587 ,  0.21707827,  0.21714922])