如何从lmfit minimizer获取拟合参数的不确定性？

使用lmfit拟合洛伦兹曲线，因参数边界限制无法直接使用scipy.optimize，拟合结果正常，但始终无法获取参数误差：res.covar返回None，res.errorbars为False，调用res.uvars触发AttributeError，res.params['parameter_name'].stderr也为None。已安装numdifftools，环境为lmfit 1.3.2、numdifftools 0.9.41、Python 3.11.7，拟合过程无警告。

拟合代码 #

def to_minimize(params, x, y, uncertainties=None):
    x0 = params["x0"]
    ampl0 = params["ampl0"]
    ampl = params["ampl"]
    sig = params["sig"]

    model = lorentzian(x, x0, ampl0, ampl, sig)
    if uncertainties is not None:
        return (y- model) / uncertainties
    else:
        return y - model


params = lmfit.Parameters()
params.add("x0", value=guess[0], min=guess[0] - 1000, max=guess[0] + 1000)
params.add("ampl0", value=guess[1], min=0)
params.add("ampl", value=guess[2], min=0)
params.add("sig", value=guess[3], min=0)

out = lmfit.minimize(to_minimize, params, args=(x_data, y_data))

拟合输出 #

Fit Result
Fit Statistics
fitting method  leastsq
# function evals    156
# data points   51201
# variables 4
chi-square  9.8949e-23
reduced chi-square  1.9327e-27
Akaike info crit.   -3149414.19
Bayesian info crit. -3149378.82
Parameters
name    value   initial value   min max vary
x0  2424682.92  2424725.9731872166  2423725.97  2425725.97  True
ampl0   1.6720e-13  1.6717888361052836e-13  0.00000000  inf True
ampl    1.2557e-12  1.2557626695414776e-12  0.00000000  inf True
sig 105.300229  4.296879768371582   0.00000000  inf True

解决办法 #

1. 强制计算协方差矩阵 #

默认leastsq方法在参数接近边界时不会自动计算误差，可手动触发计算：

# 方式1：拟合时直接指定计算协方差
out = lmfit.minimize(to_minimize, params, args=(x_data, y_data), calc_covar=True)

# 方式2：拟合后手动计算
out = lmfit.minimize(to_minimize, params, args=(x_data, y_data))
out.calc_covar()

若参数接近边界导致计算失败，可尝试用数值微分工具计算：

from lmfit.minimizer import Minimizer

minimizer = Minimizer(to_minimize, params, args=(x_data, y_data))
out = minimizer.leastsq(calc_covar=True)
# 若仍失败，强制使用numdifftools计算数值协方差
out.params.add_all_covars(minimizer.calc_covar(numdifftools=True))

2. 检查参数边界限制 #

从拟合结果看，ampl0和ampl的初始值与拟合值几乎一致，且接近下限0。如果参数被边界限制，协方差矩阵无法正常计算：

• 尝试放宽参数下限（比如设为-1e-16而非严格0），允许参数在边界内小幅波动
• 检查模型是否合理，是否存在参数冗余或过度约束

3. 手动计算参数不确定性 #

如果上述方法无效，可通过以下方式手动计算：

• 参数扫描法：固定其他参数，逐个调整目标参数，找到使卡方值增加1（单参数95%置信区间对应卡方增加4）的参数范围
• 蒙特卡洛模拟：在拟合结果基础上，给参数添加正态分布噪声，重复拟合多次，统计参数分布的标准差
• 数值微分计算协方差：用numdifftools计算模型对参数的偏导数，结合残差推导协方差：

  import numdifftools as nd
  import numpy as np
  
  def model_func(params):
      return lorentzian(x_data, params[0], params[1], params[2], params[3])
  
  # 获取拟合后的参数值
  pvals = [out.params['x0'].value, out.params['ampl0'].value, 
           out.params['ampl'].value, out.params['sig'].value]
  # 计算雅可比矩阵
  jac = nd.Jacobian(model_func)(pvals)
  # 计算残差
  resids = y_data - model_func(pvals)
  # 推导协方差矩阵
  covar = np.linalg.inv(jac.T @ jac) * (np.sum(resids**2)/(len(y_data)-4))
  # 参数标准差为协方差矩阵对角线的平方根
  stderrs = np.sqrt(np.diag(covar))
  

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