scipy.optimize 的参数拟合错误

Question

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

上面是 returns 的示例：

      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

这实际上取决于您所说的“错误”是什么意思。您的问题没有通用的答案，因为这取决于您要拟合什么以及您所做的假设。

最简单的情况也是最常见的情况之一：当您要最小化的函数是负对数似然时。在这种情况下，拟合返回的 hessian 矩阵的逆 (hess_inv) 是描述高斯近似的协方差矩阵，最大 likelihood.The 参数误差是协方差矩阵对角线元素的平方根.

请注意，如果您要拟合不同类型的函数或做出不同的假设，则这不适用。

Answer 2

解决此常见问题的一种方法是在使用 minimize 和 'L-BFGS-B' 之后使用 scipy.optimize.leastsq，从 'L-BFGS-B' 找到的解决方案开始。也就是说，leastsq 将（通常）包括和估计 1-sigma 误差以及解。

当然，该方法有几个假设，包括可以使用 leastsq 并且可能适合解决问题。从实际的角度来看，这需要 objective 函数 return 一个残值数组，其元素至少与变量一样多，而不是成本函数。

您可能会发现 lmfit (https://lmfit.github.io/lmfit-py/) 在这里很有用：它支持 'L-BFGS-B' 和 'leastsq' 并为这些和其他最小化方法提供统一的包装，因此您可以对这两种方法使用相同的 objective 函数（并指定如何将残差数组转换为成本函数）。此外，参数边界可用于这两种方法。这使得首先使用 'L-BFGS-B' 然后使用 'leastsq' 进行拟合非常容易，使用 'L-BFGS-B' 中的值作为起始值。

Lmfit 还提供了更详细地更明确地探索参数值置信限度的方法，以防您怀疑 leastsq 使用的简单但快速的方法可能不够。

Answer 3

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 作为卡方值。所以这粗略地表明你会期望

当然，您正在应用这样的最小化例程这一事实假设您的函数表现良好，在某种意义上它相当平滑并且找到的最佳值非常接近near参数 xⁱ:

的二次函数的最优值

其中Δxⁱ是参数x^{i的求值差}及其最优值，H_ij为Hessian矩阵。一点（令人惊讶的非平凡的）线性代数可以让你得到一个非常标准的结果来估计任何数量的不确定性 X 这是你的参数 x[=71 的函数=]i:

这让我们写

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

最后，为了完全明确，假设您已将优化结果存储在名为 res 的变量中。反 Hessian 可用作 res.hess_inv，它是一个函数，它接受一个向量和 returns 反 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，您显然需要在此处进行更改。

scipy.optimize 的参数拟合错误

Errors to fit parameters of scipy.optimize

numpy

scipy

curve-fitting