基于scipy.optimize的指数衰减曲线拟合优化问题

我正在尝试使用SciPy中的curve_fit函数进行指数衰减曲线拟合，调整模型初始值会改变拟合质量，但始终无法找到适配数据的最优拟合效果。

我的问题：

• 如何提升拟合质量？
• 选择模型初始值的最佳方法是什么？

原始数据 #

y = [ 338.52656636  337.43934446  348.25434126  308.42768639  279.24436171
  269.85992004  279.24436171  249.25992615  239.53215125  219.96215705
  220.41993469  220.30549028  220.30549028  195.07049776  180.364391
  171.20883816  180.24994659  180.13550218  180.47883541  209.89104892
  220.19104587  180.02105777  595.45426801  324.50712607  150.60884426
  170.97994934  171.20883816  170.75106052  170.75106052  159.76439711
  140.88106937  150.37995544  140.88106937 1620.70451979  140.42329173
  150.37995544  140.53773614  284.68047121 1146.84743797  170.97994934
  150.60884426  145.74495682  141.10995819  121.53996399  121.19663076
  131.38218329  170.40772729  140.42329173  140.82384716  145.5732902
  140.30884732  121.53996399  700.39979247 2783.74584185  131.26773888
  140.76662496  140.53773614  121.76885281  126.23218482  130.69551683]

拟合代码 #

from numpy import arange, array
import numpy as np
from pandas import read_csv
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt


def expDecay(t, Amax, tau):
    return Amax/tau*np.exp(-t/tau)

Amax = []
Tau = []

ydata = y
x = array(range(len(y)))
xdata = x


popt, pcov = curve_fit(expDecay, x, y,
                           p0=(10000, 5),
                           bounds=([0., 2.], [10000., 30]),)
Amax.append(popt[0])
Tau.append(popt[1])


plt.plot(xdata, expDecay(xdata, *popt), 'k-', label='Pred.');
plt.plot(ydata);
plt.ylim([0, 500])   
plt.show()

一、提升拟合质量的方法 #

1. 处理数据异常值 #

数据中存在多个突发高值点（如1620.70、1146.85、2783.75），这类值完全不符合指数衰减的趋势，会严重干扰拟合结果：

• 先去掉plt.ylim([0, 500])，可视化全范围数据定位异常点；
• 根据业务逻辑选择：直接剔除异常点，或通过curve_fit的sigma参数给异常点分配大权重值，降低其对拟合的影响。

2. 修正模型形式 #

当前模型假设t→∞时y趋近于0，但你的数据后期平稳在120-140区间，需加入稳态基线项：

def expDecay(t, Amax, tau, y0):
    # y0为稳态基线值，匹配数据后期平稳趋势
    return y0 + (Amax - y0)/tau * np.exp(-t/tau)

3. 调整参数约束边界 #

当前设置的tau上限为30，从数据趋势看衰减到稳态的时间更长，过窄的边界会限制拟合结果。可以先移除bounds参数观察拟合结果，再根据实际情况添加合理约束。

4. 加权拟合（保留异常点时） #

给异常点分配更低的权重，示例代码：

sigma = np.ones_like(y)
# 假设异常点索引为33、38、53，设置大sigma降低权重
sigma[[33,38,53]] = 10
popt, pcov = curve_fit(expDecay, x, y, p0=(350, 15, 130), sigma=sigma)

二、选择模型初始值的最佳方法 #

1. 从数据趋势手动估算 #

• Amax：取数据初始峰值，前几个点在330-350区间，初始值设为350即可，无需设为远高于实际的10000；
• tau：衰减时间常数，估算数据从初始值降到稳态值63%左右的时间，大概10-15个x轴单位，初始值设为15；
• y0：取数据后期平稳区间的平均值，约130。

2. 对数转换做线性拟合估算 #

把非线性指数模型转为线性模型，快速获取初始参数：

# 取前15个下降趋势明显的点
x_sub = x[:15]
y_sub = y[:15]
# 转换为ln(y) = ln(Amax/tau) - t/tau，拟合线性模型
log_y = np.log(y_sub)
a, b = np.polyfit(x_sub, log_y, 1)
tau_init = -1/b
Amax_init = np.exp(a) * tau_init

3. 网格搜索最优初始值 #

手动估算困难时，遍历多组初始值组合，选拟合误差最小的作为初始点：

from itertools import product

# 定义参数候选范围
Amax_candidates = [300, 350, 400]
tau_candidates = [10, 15, 20]
y0_candidates = [120, 130, 140]

best_err = float('inf')
best_p0 = None

for p0 in product(Amax_candidates, tau_candidates, y0_candidates):
    try:
        popt_temp, _ = curve_fit(expDecay, x, y, p0=p0)
        y_pred = expDecay(x, *popt_temp)
        err = np.sum((y - y_pred)**2)
        if err < best_err:
            best_err = err
            best_p0 = p0
    except RuntimeError:
        continue

# 使用最优初始值拟合
popt, pcov = curve_fit(expDecay, x, y, p0=best_p0)

内容的提问来源于stack exchange，提问作者Jason