利用scipy.OptimizeResult评估非线性回归中变量的相关性

用scipy.minimize结果估计参数相关性 #

完全可以通过minimize返回的OptimizeResult对象量化参数间的相关性，核心是利用参数协方差矩阵推导相关系数：

具体步骤 #

1. 获取逆海森矩阵：若使用BFGS、L-BFGS-B等带二阶信息的优化算法，结果的hess_inv属性存储了近似的逆海森矩阵。
2. 计算残差方差：用优化得到的残差平方和估计噪声方差：

   import numpy as np
   n = len(t)  # 观测样本数量
   p = len(res.x)  # 待估参数个数
   residual_var = res.fun * 2 / (n - p)  # res.fun是残差平方和的一半，需乘2还原
   

3. 推导协方差矩阵：参数协方差矩阵 = 逆海森矩阵 × 残差方差：

   cov_matrix = res.hess_inv * residual_var
   

4. 计算相关系数矩阵：通过协方差矩阵归一化得到相关系数：

   stds = np.sqrt(np.diag(cov_matrix))
   corr_matrix = cov_matrix / np.outer(stds, stds)
   

   矩阵中corr_matrix[i,j]就是第i和第j个参数的相关系数，绝对值越接近1，参数间相关性越强。

无导数算法的处理 #

如果用的是Nelder-Mead这类无导数优化算法，hess_inv会缺失，可通过数值方法计算雅可比矩阵，再用高斯-牛顿法近似协方差矩阵：

from scipy.optimize import approx_fprime

# 数值计算雅可比矩阵
jac = approx_fprime(res.x, my_func, epsilon=1e-8)
# 近似协方差矩阵
cov_matrix = residual_var * np.linalg.inv(jac.T @ jac)

关于curve_fit的args参数误解 #

curve_fit其实支持args关键字参数，正确用法是把固定常数a,b,c放在模型函数的参数列表末尾，再通过args元组传入：

def fit_model(t, x1, x2, x3, a, b, c):
    # 你的非线性模型表达式，示例：y = x1*np.exp(-x2*t) + x3*t + a*b + c
    return ...

# 调用curve_fit，传入待估参数初始值和固定常数
popt, pcov = scipy.optimize.curve_fit(
    fit_model, t, y, 
    p0=[2000.0, 0.0, 0.1], 
    args=(a, b, c)
)

这样curve_fit会自动区分待估参数（p0对应的x1,x2,x3）和固定常数（args传入的a,b,c）。

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