如何通过scipy.minimize获取MLE标准误？求可返回Hessian的优化方法

嘿，这个问题我刚好在做约束MLE估计的时候碰到过！你说得完全对——能不能拿到Hessian矩阵（进而算出标准误），全看你选的优化方法。下面给你掰扯清楚：

1. 直接能拿到Hessian/其逆的约束优化方法 #

约束优化里，最省心的是这俩方法：

• method='trust-constr': 这是scipy专门为约束优化设计的“优等生”，收敛后会在结果对象的hess_inv属性里给出Hessian的逆矩阵（或者可靠的近似逆）。如果需要精确Hessian，还能通过设置hessp参数（提供Hessian-vector乘积函数）让它计算更准确的结果。
• method='Newton-CG': 这个基于牛顿法的方法，会在优化过程中迭代计算Hessian，不过它返回的hess_inv是用共轭梯度法得到的近似逆，不是原始Hessian，但用来算标准误足够用了。

2. 其他约束方法（比如SLSQP）怎么凑Hessian？ #

如果你的问题更适合用SLSQP（这是约束优化里非常常用的方法，但默认不返回Hessian），也有两种靠谱的方式：

手动调用Hessian计算函数 #

在优化得到的最优参数x_opt处，直接用scipy.optimize.hessian算：

from scipy.optimize import minimize, hessian

# 假设你的负对数似然函数是neg_log_likelihood，约束是constraints
opt_result = minimize(neg_log_likelihood, x0, method='SLSQP', constraints=constraints)
x_best = opt_result.x

# 在最优参数处计算Hessian
hess_matrix = hessian(neg_log_likelihood, x_best)

基于梯度做数值近似 #

如果你的目标函数已经提供了梯度（也就是传了jac参数），可以用有限差分法对梯度再求一次差分，近似得到Hessian。比如用scipy.optimize.approx_fprime来算梯度的差分，代码大概是这样：

from scipy.optimize import approx_fprime
import numpy as np

# 假设grad_func是你的梯度函数
def hessian_approx(x):
    return np.array([approx_fprime(x, lambda xi: grad_func(xi)[i], epsilon=1e-6) for i in range(len(x))])

hess_matrix = hessian_approx(x_best)

最后一步：从Hessian算标准误 #

记住哦——如果你的目标函数是负对数似然（MLE里通常这么设，因为要转成最小化问题），Hessian是正定的，直接取逆然后开根号对角线元素就行：

import numpy as np

hess_inv = np.linalg.inv(hess_matrix)
# 标准误就是逆矩阵对角线的平方根
std_errors = np.sqrt(np.diag(hess_inv))

要是你用的是对数似然（最大化问题，转最小化加了负号的话其实和上面一样），那要先对Hessian取负再求逆，再算标准误。

小提醒 #

• 约束优化优先选trust-constr，它对约束的处理更稳健，返回的Hessian近似也更可靠；
• 手动算Hessian的时候，最好检查一下Hessian是不是正定的——如果不是，要么是优化没收敛到全局最优，要么是你的模型设定有问题；
• 要是优化方法直接返回了hess_inv（比如trust-constr），直接用它就行，不用再手动求逆啦。

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