使用curve_fit拟合负指数模型时L_infty参数不确定性过高的优化方法

问题背景 #

使用负指数模型twentyconc_model拟合5组数据，数据预期符合L_infty*(1 - exp(-k*x))的饱和增长关系，但拟合得到的饱和值L_infty预测不确定性远大于参数k，误差棒显示其波动极大，需要优化拟合效果。

核心原因分析 #

1. 参数相关性：模型中L_infty与k存在较强的共线性，当数据未充分饱和时，拟合算法难以同时准确推断两个参数，导致L_infty的不确定性被放大。
2. 初始值不合理：k的初始值设为0会拖慢拟合收敛，若数据未完全饱和，仅用末尾10个点的均值作为L_infty初始值也不够准确。
3. 无参数约束：未对L_infty（需大于数据最大值）、k（需为正数）设置边界，拟合过程中参数可能进入不合理区间。
4. 数据范围不足：若数据的x（时间）范围不够大，未接近饱和阶段，模型无法准确推断L_infty的真实值。

具体优化措施 #

1. 添加参数边界约束 #

利用curve_fit的bounds参数，强制L_infty大于数据中的最大值（饱和值不可能低于数据观测值），k为正数，避免参数进入不合理区间：

# 设置边界：k>0，L_infty不小于数据最大值的1.01倍，上限可根据实际情况调整
bounds = ([0, np.max(y_data)*1.01], [np.inf, np.max(y_data)*1.5])

2. 优化初始值设置 #

• 避免k初始值为0，通过数据初始斜率估算合理的正数初始值：

  # 用前两个点的斜率估算k的初始值
  Linfty_guess = np.mean(y_data[-5:]) if len(y_data)>=5 else np.max(y_data)
  initial_slope = (y_data[1] - y_data[0])/(x_data[1]-x_data[0])
  k_guess = initial_slope / Linfty_guess if Linfty_guess !=0 else 0.01
  

• 若数据未饱和，用末尾点的线性外推值替代简单均值作为L_infty初始值。

3. 重构模型或联合拟合（基于物理规律） #

如果多组数据（如不同温度下的实验）存在已知物理关联，比如L_infty随温度线性变化，可将多组数据合并拟合，共享参数关系，减少自由度：

# 示例：假设温度与L_infty线性相关，Linfty = a*T + b
def combined_model(x, k, a, b, T):
    Linfty = a*T + b
    return Linfty*(1 - np.exp(-k*x))

# 合并所有数据进行拟合
all_x = np.concatenate([twentyconc_data[i][:,0] for i in range(5)])
all_y = np.concatenate([twentyconc_data[i][:,1] for i in range(5)])
all_T = np.concatenate([np.full(len(twentyconc_data[i]), T_list[i]) for i in range(5)])  # T_list为各组温度

popt_combined, pcov_combined = curve_fit(combined_model, all_x, all_y, p0=[0.01, 0.1, 1.0], bounds=([0, -np.inf, -np.inf], [np.inf, np.inf, np.inf]))

4. 数据预处理与补充 #

• 对噪声较大的数据做平滑处理（如滑动平均），减少噪声对拟合的干扰；
• 若数据未接近饱和，延长实验时间获取更接近饱和的观测值，这是降低L_infty不确定性最直接的方法。

调整后的完整拟合代码 #

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

def twentyconc_model(x, k, Linfty):
    return Linfty * (1 - np.exp(-k * x))  # 等价于原模型，写法更清晰

popt20 = [] 
pcov20 = []
label20 = ["组1", "组2", "组3", "组4", "组5"]  # 替换为实际标签

for i in range(5):
    x_data = twentyconc_data[i][:, 0]
    y_data = twentyconc_data[i][:, 1]
    
    # 估算初始值
    Linfty_guess = np.mean(y_data[-5:]) if len(y_data)>=5 else np.max(y_data)
    k_guess = 0.01  # 默认初始值
    if len(x_data) >= 2:
        initial_slope = (y_data[1] - y_data[0]) / (x_data[1] - x_data[0])
        if Linfty_guess != 0:
            k_guess = initial_slope / Linfty_guess
        # 确保k初始值为正
        k_guess = max(k_guess, 0.001)
    
    # 设置参数边界
    y_max = np.max(y_data)
    bounds = ([0, y_max * 1.01], [np.inf, y_max * 1.5])
    
    # 拟合
    popt, pcov = curve_fit(
        twentyconc_model,
        x_data,
        y_data,
        p0=[k_guess, Linfty_guess],
        maxfev=10000,
        bounds=bounds
    )
    popt20.append(popt)
    pcov20.append(pcov)

# 输出结果
print("各组拟合参数(k, Linfty):", popt20)
print("各组L_infty的标准差(不确定性):", [np.sqrt(pcov[1,1]) for pcov in pcov20])

# 绘制拟合曲线与原始数据
plt.figure(figsize=(10, 6))
for i in range(5):
    x_data = twentyconc_data[i][:, 0]
    y_data = twentyconc_data[i][:, 1]
    k_fit, Linfty_fit = popt20[i]
    
    # 绘制原始数据点
    plt.scatter(x_data, y_data, s=20, label=f"{label20[i]} 原始数据")
    # 绘制拟合曲线
    x_model = np.linspace(min(x_data), max(x_data), 100)
    y_model = twentyconc_model(x_model, k_fit, Linfty_fit)
    plt.plot(x_model, y_model, linewidth=2, label=f"{label20[i]} 拟合曲线")

plt.xlabel("时间")
plt.ylabel("电导")
plt.legend()
plt.grid(alpha=0.3)
plt.show()

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