使用scipy.optimize.curve_fit拟合高斯包络余弦曲线失败排查

问题描述 #

尝试用scipy.optimize.curve_fit拟合带衰减高斯包络的余弦函数数据时，出现参数优化不充分的问题：干涉图中心水平偏移量（约-0.1）的拟合结果始终紧贴近初始猜测值，而非手动验证的最优值，即便初始猜测已经接近真实值。手动调整水平偏移量并计算决定系数后，确认curve_fit给出的结果并非最优解。拟合结果显示，拟合曲线与原始数据在中心位置存在明显偏移。

拟合失败的核心原因 #

• 参数强耦合：拟合函数FF_function中，水平偏移量xoff与余弦项周期l存在强相关性。当xoff变化时，微调l即可让曲线在数据上表现出近似的拟合效果，导致优化算法陷入局部极小值，无法收敛到全局最优的xoff。
• 参数数量级差异过大：各参数数量级跨度明显（如A为100+，l为0.04+，xoff为-0.2），curve_fit默认的Levenberg-Marquardt算法对参数缩放敏感，易导致部分参数更新步长不合理，难以找到最优解。
• 基线设定不合理：函数中基线固定为A，但实际数据的基线未必严格等于振幅A，额外的拟合误差会干扰算法对xoff等参数的优化。

优化建议 #

1. 重构拟合函数，解耦参数
   引入独立的基线参数B，避免基线与振幅绑定，同时显式建模相位，减少xoff与l的耦合：

   def FF_function(x, l, s, A, B, xoff):
       d = x - xoff
       return A * np.cos(2 * np.pi / l * d) * np.exp(-0.5 * (d / s)**2) + B
   

2. 更换鲁棒性更强的优化算法
   使用method='trf'（Trust Region Reflective）算法，它对参数缩放的鲁棒性优于默认的Levenberg-Marquardt算法：

   popt, cov = curve_fit(FF_function, data[:,0], data[:,1], p0=guess, method='trf')
   

3. 加入参数边界约束
   根据数据先验知识设置参数上下界，避免算法在无意义的参数空间搜索：

   # 示例边界，需根据实际数据调整
   bounds = (
       [0.04, 0.2, 100, 0, -0.3],  # 参数下界
       [0.05, 0.3, 200, 500, 0.1]  # 参数上界
   )
   

4. 优化初始猜测值

   • 用数据均值作为基线B的初始值
   • 通过数据峰值/谷值位置手动估算更准确的xoff初始值，减少算法优化负担

修改后的完整代码 #

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

id_ = '0' 
filename = '250516_FF_01_coincidences_'+id_+'mm.txt'

# 重构拟合函数，加入独立基线参数B
def FF_function(x, l, s, A, B, xoff):
    d = x - xoff 
    return A * np.cos(2 * np.pi / l * d) * np.exp(-0.5 * (d / s)**2) + B

data = pd.read_csv(filename, sep=' ', header=None).to_numpy()

fig, ax = plt.subplots(figsize=(10,5))
ax.plot(data[:,0], data[:,1], alpha=0.3, label='原始数据')

# 更新初始猜测：加入基线B的初始值（数据均值），优化xoff初始值
guess = [4.4067e-02, 2.6397e-01, 1.6370e+02, np.mean(data[:,1]), -0.1]
xx = np.linspace(min(data[:,0]), max(data[:,0]), 1000)

# 配置参数边界并使用trf算法
bounds = (
    [0.04, 0.2, 100, 0, -0.3],
    [0.05, 0.3, 200, 500, 0.1]
)
popt, cov = curve_fit(FF_function, data[:,0], data[:,1], p0=guess, method='trf', bounds=bounds)
ax.plot(xx, FF_function(xx, *popt), color='C0', alpha=0.8, label='拟合曲线')

# 计算决定系数
SSE = np.sum(np.square(data[:,1] - FF_function(data[:,0], *popt)))
SST = np.sum(np.square(data[:,1] - np.mean(data[:,1])))
r2 = 1 - SSE/SST
print(f"决定系数R²: {r2:.4f}")
print(f"拟合参数: {popt}")

ax.set_xlabel('Displacement (mm)')
ax.set_ylabel('Coincidences')
ax.legend()
plt.show()

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