Python中加权最小二乘拟合异常问题求助

问题背景 #

需要对给定数据进行加权最小二乘拟合，要求：

• x=0和x>0.2的数据点权重较高
• 0<x≤0.2的点权重较低但保留一定权重
  当前拟合效果差（大x值偏差明显），且出现divide by zero encountered in power运行警告，期望拟合曲线趋近于参考图（曲线在x=0处为0，x>0.2时平稳趋近于常数，中间区间贴合数据趋势）。

原始数据 #

import numpy as np

filtered_x = np.array([0.        , 0.03807991, 0.05077321, 0.06346652, 0.07615982,
       0.08885313, 0.10154643, 0.11423973, 0.12693304, 0.13962634,
       0.15231964, 0.16501295, 0.17770625, 0.19039955, 0.20309286,
       0.21578616, 0.22847947, 0.24117277, 0.25386607, 0.26655938,
       0.27925268, 0.29194598, 0.30463929, 0.31733259, 0.33002589,
       0.3427192 ])

filtered_y = np.array([0.        , 1.53989397, 2.04460628, 4.18043213, 2.97621482,
       2.82642339, 2.98335023, 2.98964836, 2.12218901, 1.42801972,
       1.25930683, 0.71644077, 0.48220866, 0.21165985, 0.24756609,
       0.21123179, 0.57344999, 0.49362762, 0.20282767, 0.50321186,
       0.50347165, 0.74259408, 0.48493783, 0.81785588, 0.54543666,
       0.53218838])

现有代码问题分析 #

1. 拟合函数的零值处理缺陷：当x=0且拟合参数b为负数时，0**b会触发除以零的计算错误。
2. 初始猜测值不合理：原初始值(0.4, 0.4, 1)不符合数据趋势（参考图显示曲线最终平稳在~0.5，且b应为负数）。
3. 拟合收敛性不足：虽然设置了maxfev=5000，但初始值偏差过大可能导致拟合未收敛到最优解。

修正后的完整代码 #

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

# 修正后的拟合函数：处理x=0的特殊情况，避免除以零错误
def fpow(x, a, b, c):
    x_arr = np.asarray(x)
    # x=0时直接返回c，保证该点计算合法且贴合数据
    return np.where(x_arr == 0, c, a * (x_arr ** b) + c)

# 权重函数（逻辑正确：sigma越小，点权重越高）
def custom_weights(x):
    weights = np.ones_like(x)
    # 高权重点（误差小）：x=0和x>0.2
    weights[x == 0] = 0.1
    weights[x > 0.2] = 0.1
    # 低权重点（误差大）：0<x≤0.2
    weights[(0 < x) & (x <= 0.2)] = 0.5
    return weights

# 调整初始猜测值：贴合数据趋势，b设为负数，c接近x>0.2的y均值
pars0 = (10, -1, 0.5)

# 执行加权拟合
popt, pcov = curve_fit(
    fpow, 
    filtered_x, 
    filtered_y, 
    absolute_sigma=True, 
    p0=pars0, 
    sigma=custom_weights(filtered_x), 
    maxfev=10000  # 增加迭代次数保证收敛
)

# 提取最优参数
a_opt, b_opt, c_opt = popt
print(f"最优参数：a={a_opt:.3f}, b={b_opt:.3f}, c={c_opt:.3f}")

# 绘图展示
plt.errorbar(filtered_x, filtered_y, yerr=0, fmt=".", color='black', label='data', zorder=1, markersize=5)

# 生成拟合曲线的x区间
x_interval = np.linspace(0, max(filtered_x), 1000)
y_fit = fpow(x_interval, *popt)
plt.plot(x_interval, y_fit, color="red", label="Weighted fit", zorder=2, linewidth=3)

plt.grid(True)
plt.ylabel("U [m/s]")
plt.xlabel("Wavenumber [rad/m]")
plt.title("LS Weighted Fit of Current")
plt.legend()
plt.show()

关键修正说明 #

1. 拟合函数修复：通过np.where处理x=0的情况，彻底消除除以零警告，同时保证x=0处的拟合值符合预期。
2. 初始值优化：将b设为负数（符合曲线下降趋势），c设为0.5（接近x>0.2时的y值均值），大幅提升拟合收敛效率和准确性。
3. 迭代次数增加：将maxfev提升至10000，确保拟合算法有足够次数收敛到最优解。

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