最小化拟合解的误差为何过大？理论模型拟合求助

问题概述 #

我有若干大型数据集，正在构建理论模型进行拟合。已找到全局最小值，chi² landscape清晰显示~400附近存在明确最小值，但计算目标参数误差时结果远大于预期。模型包含1个核心目标可调参数，以及3个用于全局修正模型与实验数据匹配的缩放因子。

拟合实现代码 #

import numpy as np
from scipy.optimize import minimize

def fun(k,io,exp_data,exp_error,model1,model2):
    F=(np.sqrt((8*k[0]*io)+k[0]**2)-k[0])/4
    C=(-np.sqrt(k[0])*np.sqrt((8*io)+k[0])+k[0]+(4*io))/4
    fraction_F=np.tile(((np.array([F/io],dtype=float)).reshape((len(model1),1))),(len(model1[0]))) # 单值应用于所有数据集
    fraction_C=np.tile(((np.array([C/io],dtype=float)).reshape((len(model1),2))),(len(model2[0])))
    comb_model=(F*model1)+(C/2*model2)
    resize_scaling_factor=np.tile(((np.array(k[1:])).reshape((len(comb_model),1))),(len(comb_model[0]))) # 全局缩放因子应用于数据集
    return ((np.sum(((exp_data-(resize_scaling_factor*comb_model))/exp_error)**2))/comb_model.size) # 归一化chi²

minimize(fun,args=(io,exp_data,exp_error,model1,model2),bounds=((0,np.inf),)*4,method='Nelder-Mead',x0=[1000]+[1e-9]*3)

拟合结果与矛盾点 #

• 拟合能顺利收敛：目标参数解在500-1000区间，缩放因子为1e-8至1e-11量级
• chi² landscape显示最小值明确，按经验目标参数误差应在±50左右，但实际计算得到：

  428+/-593 # 第一个值是极小化解，第二个是逆Hessian对角线开方结果
  # 缩放因子误差表现正常，相对解的比例很小
  7.91e-10+/-4.5e-11
  1.72e-9+/-5.34e-11
  3.52e-9+/-6.3e-11 
  

• 归一化chi²接近1，符合拟合质量预期，但目标参数误差异常偏大

误差计算代码 #

from numpy.linalg import inv
from statsmodels.tools import approx_hess2
solution=minimize(fun,args=(io,exp_data,exp_error,model1,model2),bounds=((0,np.inf),)*4,method='Nelder-Mead',x0=[1000]+[1e-9]*3)
error=np.sqrt(np.diag(inv(approx_hess2(solution.x,fun,args=(io,exp_data,exp_error,model1,model2)))))

可能的原因与解决方向 #

1. 参数尺度不匹配导致Hessian数值不稳定 #

目标参数（500-1000）与缩放因子（1e-8至1e-11）量级相差11-14个数量级，会导致Hessian矩阵条件数极差，求逆时数值不稳定，放大目标参数的误差估计。

• 解决方法：对目标参数做尺度归一化（例如除以1000，转化为0.5-1区间），拟合后再还原参数与误差：

  # 修改目标参数的缩放逻辑
  def fun_scaled(k_scaled,io,exp_data,exp_error,model1,model2):
      k0 = k_scaled[0] * 1000  # 还原为原尺度参数
      k = [k0] + list(k_scaled[1:])
      # 原拟合逻辑不变
      F=(np.sqrt((8*k[0]*io)+k[0]**2)-k[0])/4
      C=(-np.sqrt(k[0])*np.sqrt((8*io)+k[0])+k[0]+(4*io))/4
      fraction_F=np.tile(((np.array([F/io],dtype=float)).reshape((len(model1),1))),(len(model1[0])))
      fraction_C=np.tile(((np.array([C/io],dtype=float)).reshape((len(model1),2))),(len(model2[0])))
      comb_model=(F*model1)+(C/2*model2)
      resize_scaling_factor=np.tile(((np.array(k[1:])).reshape((len(comb_model),1))),(len(comb_model[0])))
      return ((np.sum(((exp_data-(resize_scaling_factor*comb_model))/exp_error)**2))/comb_model.size)
  
  # 初始值对应缩放后的尺度
  solution_scaled=minimize(fun_scaled,args=(io,exp_data,exp_error,model1,model2),bounds=((0,np.inf),)*4,method='Nelder-Mead',x0=[1.0]+[1e-9]*3)
  # 还原参数
  solution_x = [solution_scaled.x[0]*1000] + list(solution_scaled.x[1:])
  # 计算Hessian并还原误差
  hess_scaled = approx_hess2(solution_scaled.x,fun_scaled,args=(io,exp_data,exp_error,model1,model2))
  error_scaled = np.sqrt(np.diag(inv(hess_scaled)))
  error = [error_scaled[0]*1000] + list(error_scaled[1:])
  

2. 无导数优化方法的Hessian近似精度不足 #

Nelder-Mead是无导数优化方法，最优解附近的梯度信息精度有限，叠加数值Hessian的近似误差，在参数尺度差异大的场景下会被放大。

• 解决方法：改用带导数的优化方法（如L-BFGS-B），若能手动推导目标函数梯度，精度会进一步提升；或使用scipy.optimize.least_squares，该工具专门针对最小二乘问题，能更可靠地估计参数协方差。

3. 归一化chi²的尺度影响 #

目标函数返回的是归一化chi²（除以数据点数量），会改变Hessian矩阵的尺度，导致逆Hessian的误差估计需要调整。标准最小二乘中，参数协方差矩阵是inv(Hessian)乘以归一化chi²，直接使用归一化目标函数可能导致误差尺度错误。

• 解决方法：尝试将目标函数改为未归一化的chi²（移除/comb_model.size），重新计算Hessian与误差，观察结果是否符合预期。

4. 数值Hessian的步长不匹配 #

approx_hess2的默认步长不适合量级差异极大的参数，会导致Hessian近似误差过大。

• 解决方法：手动指定对应参数的步长，对目标参数用较大步长（如10），对缩放因子用极小步长（如1e-12）：

  epsilon = [10, 1e-12, 1e-12, 1e-12]  # 对应每个参数的计算步长
  hess = approx_hess2(solution.x,fun,args=(io,exp_data,exp_error,model1,model2),epsilon=epsilon)
  error=np.sqrt(np.diag(inv(hess)))
  

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