高斯牛顿法拟合温度-时间数据参数失败求助

核心错误分析 #

拟合效果极差主要由三个关键问题导致：

• 雅可比矩阵（偏导数）计算错误：对参数b和c的偏导数完全搞反，且对b的偏导缺少1/c的因子，导致Gauss-Newton法的方向更新完全偏离最优路径。
• 迭代更新符号错误：错误使用x = x - s进行参数更新，违背Gauss-Newton法的参数更新逻辑。
• 初始参数选择不合理：d的初始值设为0，与温度数据的实际均值偏差过大，迭代难以收敛到最优解。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# 加载测量温度数据
_, T_measured = np.genfromtxt("4_temperatures.txt", unpack=True)
t_measured = np.arange(T_measured.shape[0])  # 距首次测量的天数

def T(t, p):
    a, b, c, d = p
    return a * np.sin((t - b)/c) + d

def DT(t, p):
    a, b, c, d = p
    arg = (t - b)/c
    cos_arg = np.cos(arg)
    # 正确计算各参数的偏导数：∂T/∂a, ∂T/∂b, ∂T/∂c, ∂T/∂d
    return np.array([
        np.sin(arg),                  # ∂T/∂a
        -a * cos_arg / c,             # ∂T/∂b
        -a * (t - b) * cos_arg / (c**2),  # ∂T/∂c
        np.ones_like(t)               # ∂T/∂d
    ]).T

def gn(x0, T, DT, tol, max_iter=1000):
    x = np.array(x0, dtype=np.float64)
    for i in range(max_iter):
        J = DT(t_measured, x)
        r = T_measured - T(t_measured, x)
        # 最小二乘求解JΔx = r
        s, _, _, _ = np.linalg.lstsq(J, r, rcond=None)
        # 正确更新参数：x = x + Δx
        x_new = x + s
        # 收敛判断
        if np.linalg.norm(s) < tol * np.linalg.norm(x_new):
            return x_new
        x = x_new
    # 若迭代到最大次数仍未收敛，返回当前值
    return x

if __name__ == "__main__":
    # 更合理的初始参数
    d_init = np.mean(T_measured)
    a_init = (np.max(T_measured) - np.min(T_measured)) / 2
    p0 = [a_init, 100, 58, d_init]  # [a, b, c, d]
    
    p = gn(p0, T, DT, 1e-8)
    
    fit = T(t_measured, p)
    
    plt.figure()
    plt.plot(t_measured, T_measured, label='实测数据')
    plt.plot(t_measured, T(t_measured, p0), label='初始猜测')
    plt.plot(t_measured, fit, label='最优拟合')
    
    plt.legend()
    plt.xlabel('距首次测量的天数')
    plt.ylabel('温度')
    plt.savefig("4_mythenquai.png")
    
    # 输出参数
    print("初始值:")
    print("a = {:.3f}, b = {:.3f}, c = {:.3f}, d = {:.3f}".format(*p0))
    print("\n最优拟合参数:")
    print("a = {:.3f}, b = {:.3f}, c = {:.3f}, d = {:.3f}".format(*p))

关键修正说明 #

1. 修正雅可比矩阵：
   • 对b的偏导数：$\partial T/\partial b = -\frac{a}{c}\cos\left(\frac{t-b}{c}\right)$
   • 对c的偏导数：$\partial T/\partial c = -\frac{a(t-b)}{c^2}\cos\left(\frac{t-b}{c}\right)$
2. 迭代更新逻辑：将x = x - s改为x = x + s，符合Gauss-Newton法的参数更新规则。
3. 优化初始参数：
   • d初始值设为温度数据的均值，保证初始曲线贴合数据中心
   • a初始值设为数据波动范围的一半，匹配正弦曲线振幅的物理意义
   • c初始值设为58，对应一年365天的周期（周期$T=2\pi c$，$c=365/(2\pi)\approx58$）

修正后，Gauss-Newton法会快速收敛，拟合曲线将紧密贴合实测温度数据，参数结果符合温度随时间周期性变化的规律。

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