理解scipy.optimize.curve_fit：自定义拟合代码报错排查

这个错误的核心是调用l_sq函数时没有传入必需的par参数，以下是具体排查和解决步骤：

1. 检查l_sq函数定义与调用的参数匹配 #

确保函数定义的参数顺序和调用时的传参顺序完全一致：

• 若损失函数定义为：

def l_sq(par, t, y):
    # 计算残差平方和：(y - sin(ωt))^2 总和，par[0]对应ω
    return np.sum((y - np.sin(par[0] * t))**2)

调用时必须先传par，再传t和y，不能颠倒或遗漏。

• 若你把par放在最后（def l_sq(t, y, par):），则调用时必须按t, y, par的顺序传参。

2. 检查牛顿-拉夫逊迭代中的函数调用 #

在迭代过程中，计算损失值、梯度、海森矩阵时，必须给l_sq传全参数：

• 错误调用：current_loss = l_sq(t, y)（漏传par）
• 正确调用：current_loss = l_sq(par_current, t, y)（对应上述第一种函数定义）

3. 修正后的完整示例代码 #

以下是可运行的修正版代码，包含中心差分近似的牛顿-拉夫逊拟合逻辑：

import numpy as np

# 生成带噪声的模拟数据
t = np.linspace(0, 2*np.pi, 100)
true_omega = 2.5
y = np.sin(true_omega * t) + 0.1*np.random.randn(len(t))

# 最小二乘损失函数：参数par在前，数据t、y在后
def l_sq(par, t, y):
    omega = par[0]
    return np.sum((y - np.sin(omega * t))**2)

# 中心差分计算梯度
def grad_l_sq(par, t, y, h=1e-6):
    grad = np.zeros_like(par)
    for i in range(len(par)):
        par_plus = par.copy()
        par_plus[i] += h
        par_minus = par.copy()
        par_minus[i] -= h
        grad[i] = (l_sq(par_plus, t, y) - l_sq(par_minus, t, y)) / (2*h)
    return grad

# 中心差分计算海森矩阵
def hess_l_sq(par, t, y, h=1e-6):
    n = len(par)
    hess = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            par_ij_plus = par.copy()
            par_ij_plus[i] += h
            par_ij_plus[j] += h
            
            par_i_plus_j_minus = par.copy()
            par_i_plus_j_minus[i] += h
            par_i_plus_j_minus[j] -= h
            
            par_i_minus_j_plus = par.copy()
            par_i_minus_j_plus[i] -= h
            par_i_minus_j_plus[j] += h
            
            par_ij_minus = par.copy()
            par_ij_minus[i] -= h
            par_ij_minus[j] -= h
            
            hess[i,j] = (l_sq(par_ij_plus, t, y) - l_sq(par_i_plus_j_minus, t, y) - l_sq(par_i_minus_j_plus, t, y) + l_sq(par_ij_minus, t, y)) / (4*h**2)
    return hess

# 牛顿-拉夫逊迭代拟合
def newton_raphson(t, y, init_par, max_iter=100, tol=1e-6):
    par = init_par.copy()
    for i in range(max_iter):
        grad = grad_l_sq(par, t, y)
        hess = hess_l_sq(par, t, y)
        # 加小正则项避免海森矩阵奇异
        hess_reg = hess + 1e-8 * np.eye(len(par))
        delta_par = np.linalg.solve(hess_reg, -grad)
        par += delta_par
        if np.linalg.norm(delta_par) < tol:
            print(f"收敛于第 {i+1} 次迭代")
            break
    return par

# 执行拟合
init_omega = np.array([1.0])
fitted_par = newton_raphson(t, y, init_omega)
print(f"拟合得到的角频率ω: {round(fitted_par[0], 4)}")

关键修正点 #

• 所有调用l_sq的位置都传入了par、t、y三个参数，无遗漏
• 统一将待优化参数par作为损失函数的第一个参数，符合迭代逻辑的参数传递习惯

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