如何在Python中优化误差函数的峰区数据拟合效果

我正在使用包含参数a,b,mu,sigma的误差函数对数据进行拟合，但峰区（如下方高亮区域）的拟合效果不佳。请问该如何优化误差函数的拟合效果？我也可通过修改func(x, a, b, mu, sigma)来提升峰区的拟合质量。

现有拟合代码 #

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

def func(x, a, b, mu, sigma):
    y = a * scipy.special.erf((x-mu) / (sigma * np.sqrt(2))) + b
    return y

def fit_and_print_params(x, y, label, color, linestyle):
    # Initial guess for the parameters
    p0 = [1, 0, 1, 10]

    # Fit the curve
    popt, pcov = optimize.curve_fit(func, x, y, p0=p0, maxfev=8000)
    a_opt = popt[0]
    b_opt = popt[1]
    mu_opt = popt[2]
    sigma_opt = popt[3]

    print(f"Optimized values for {label}:")
    print("a =", a_opt)
    print("b =", b_opt)
    print("mu =", mu_opt)
    print("sigma =", sigma_opt)

    # Calculate R-squared
    y_fit = func(x, *popt)
    residuals = y - y_fit
    ss_residuals = np.sum(residuals**2)
    ss_total = np.sum((y - np.mean(y))**2)
    r_squared = 1 - (ss_residuals / ss_total)
    print("R-squared:", r_squared)

    plt.plot(x, func(x, *popt), label=f"Fitted Curve ({label})", color=color, linestyle=linestyle)

Folder1 = '4_1012_nodes_var_1_K_1e3_Cb_1e-3_D_1e-10'
Folder2 = '4_1012_nodes_var_5_K_1e3_Cb_1e-3_D_1e-10'
Folder3 = '4_1012_nodes_var_10_K_1e3_Cb_1e-3_D_1e-10'
Drive='D'

K=1.0e3
Cb0=0.001
DiffCoeff=1.0e-10

y1 = pd.read_csv(rf'{Drive}:\\Users\\debanikb\\OneDrive - Technion\\Research_Technion\\Python_PNM\\Surfactant A-D\\{Folder1}\\Averaged4_0.5__{K}_{Cb0}_{DiffCoeff}_501files.txt',header=None)
y1 =(1 - np.ravel(y1.to_numpy() / 1012))

y2 = pd.read_csv(rf'{Drive}:\\Users\\debanikb\\OneDrive - Technion\\Research_Technion\\Python_PNM\\Surfactant A-D\\{Folder2}\\Averaged4_0.5__{K}_{Cb0}_{DiffCoeff}_501files.txt',header=None)
y2 =(1 - np.ravel(y2.to_numpy() / 1012))

y3 = pd.read_csv(rf'{Drive}:\\Users\\debanikb\\OneDrive - Technion\\Research_Technion\\Python_PNM\\Surfactant A-D\\{Folder3}\\Averaged4_0.5__{K}_{Cb0}_{DiffCoeff}_501files.txt',header=None)
y3 =(1 - np.ravel(y3.to_numpy() / 1012))

x = pd.read_csv(rf'{Drive}:\\Users\\debanikb\\OneDrive - Technion\\Research_Technion\\Python_PNM\\Surfactant A-D\\{Folder1}\\All_values_0_1000001_10000.txt',header=None)
x = np.ravel(x.to_numpy())

plt.plot(x, y1, 'o',linewidth=0.5,label='Simulation (var 1)')
fit_and_print_params(x, y1, 'var 1', color='red', linestyle='--')

plt.plot(x, y2, 'o', label='Simulation (var 5)')
fit_and_print_params(x, y2, 'var 5', color='blue', linestyle='-.') 

plt.plot(x, y3, 'o', label='Simulation (var 10)')
fit_and_print_params(x, y3, 'var 10', color='black', linestyle=':')

plt.xlabel("Elapsed time (sec)",fontsize=15.0)
plt.ylabel("Invaded fraction",fontsize=15.0)
plt.legend(loc='lower right')
plt.xlim(0, 1e6)
plt.show()

具体优化方案 #

1. 修正初始参数猜测 #

当前初始值p0 = [1, 0, 1, 10]和实际数据偏差极大，尤其是mu（峰的中心位置）和sigma（峰的宽度）。建议从图中手动读取峰对应的x值作为mu初始值，根据峰的左右跨度估算sigma，比如如果峰在x=2e5附近，跨度约1e5，可以设置p0 = [0.5, 0.5, 2e5, 5e4]，让拟合算法更快收敛到最优解，避免陷入局部最优。

2. 引入加权拟合 #

给峰区数据点更高权重，让拟合过程更关注该区域的误差。比如判断峰区范围（如1e5 < x < 5e5），给这些点设置权重2，其他点权重1，然后在curve_fit中传入权重的倒数作为残差标准差：

weights = np.where((x > 1e5) & (x < 5e5), 2, 1)
popt, pcov = optimize.curve_fit(func, x, y, p0=p0, maxfev=8000, sigma=1/weights)

3. 修改拟合函数形式 #

单一误差函数的灵活性不足，可尝试以下改进：

• 添加峰区陡峭度参数：引入额外参数调整峰区的过渡斜率：

def func(x, a, b, mu, sigma, c):
    y = a * scipy.special.erf(c*(x-mu) / (sigma * np.sqrt(2))) + b
    return y

• 叠加双误差函数：如果峰区存在复合特征，用两个误差函数叠加拟合：

def func(x, a1, b1, mu1, sigma1, a2, b2, mu2, sigma2):
    y1 = a1 * scipy.special.erf((x-mu1)/(sigma1*np.sqrt(2))) + b1
    y2 = a2 * scipy.special.erf((x-mu2)/(sigma2*np.sqrt(2))) + b2
    return y1 + y2

• 替换为逻辑斯蒂函数：逻辑斯蒂函数的S型过渡更灵活，可能更适配你的峰区形状：

def func(x, a, b, mu, sigma):
    y = a / (1 + np.exp(-(x - mu)/sigma)) + b
    return y

4. 约束参数范围 #

通过bounds参数限制参数的物理合理范围，避免拟合出异常值：

# 示例：a在0-1，b在0-1，mu在0-1e6，sigma在1e4-1e5
bounds = ([0, 0, 0, 1e4], [1, 1, 1e6, 1e5])
popt, pcov = optimize.curve_fit(func, x, y, p0=p0, maxfev=8000, bounds=bounds)

5. 数据预处理 #

• 检查峰区是否存在异常值，剔除偏离趋势的点；
• 对x轴做对数变换（如果数据随时间呈指数变化），让峰区特征更明显：

x_log = np.log10(x[x>0])  # 避开x=0的对数无意义点
# 用x_log拟合后，再转换回原x轴绘制曲线

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