numpy - 拟合 scipy.optimize 参数的错误

Question

我将 scipy.optimize.minimize ( https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html ) 函数与 method='L-BFGS-B 一起使用。

它返回的一个例子是上面的:

      fun: 32.372210618549758
 hess_inv: <6x6 LbfgsInvHessProduct with dtype=float64>
     jac: array([ -2.14583906e-04,   4.09272616e-04,  -2.55795385e-05,
         3.76587650e-05,   1.49213975e-04,  -8.38440428e-05])
  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 420
      nit: 51
   status: 0
  success: True
        x: array([ 0.75739412, -0.0927572 ,  0.11986434,  1.19911266,  0.27866406,
       -0.03825225])

x 值正确包含拟合参数。如何计算与这些参数相关的错误？

Answer 1

TL;DR: 实际上，您可以为最小化例程找到参数最佳值的精确度设置上限。请参阅此答案末尾的片段，其中显示了如何直接执行此操作，而无需调用其他最小化例程。

这种方法的 documentation 说

The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol.

粗略地说，当您要最小化的函数 f 的值最小化到最优值的 ftol 以内时，最小化停止。 (如果 f 大于 1，这是一个相对错误，否则是绝对错误；为简单起见，我假设它是一个绝对错误。)在更标准的语言中，您可能会将函数 f 视为卡方值。所以这大致表明你会期望

当然，您正在应用这样的最小化例程这一事实假设您的函数表现良好，因为它相当平滑，并且找到的最佳值通过参数 xi 的二次函数在最佳值附近很好地逼近:

其中 Δxi 是参数 xi 的发现值与其最优值之间的差值，Hij 是 Hessian 矩阵。一点点(令人惊讶的是非平凡的)线性代数让你得到一个非常标准的结果，用于估计任何数量 X 的不确定性，这是你的参数 xi 的函数:

这让我们写

这是一般最有用的公式，但对于这里的具体问题，我们只有 X = xi，所以这简化为

最后，为了完全明确，假设您已将优化结果存储在名为 res 的变量中。逆 Hessian 可用作 res.hess_inv ，这是一个函数，它接受一个向量并返回逆 Hessian 与该向量的乘积。因此，例如，我们可以使用如下片段显示优化参数以及不确定性估计:

ftol = 2.220446049250313e-09
tmp_i = np.zeros(len(res.x))
for i in range(len(res.x)):
    tmp_i[i] = 1.0
    hess_inv_i = res.hess_inv(tmp_i)[i]
    uncertainty_i = np.sqrt(max(1, abs(res.fun)) * ftol * hess_inv_i)
    tmp_i[i] = 0.0
    print('x^{0} = {1:12.4e} ± {2:.1e}'.format(i, res.x[i], uncertainty_i))

请注意，我已经合并了文档中的 max 行为，假设 f^k 和 f^{k+1} 基本上与最终输出值 res.fun 相同，这确实应该是一个很好的近似值。此外，对于小问题，您可以使用 np.diag(res.hess_inv.todense()) 来获得完整的逆并一次性提取对角线。但是对于大量变量，我发现这是一个慢得多的选择。最后，我添加了 ftol 的默认值，但是如果您在 minimize 的参数中更改它，您显然需要在此处更改它。