适合 scipy.optimize 参数的错误答案

【问题标题】：Errors to fit parameters of scipy.optimize适合 scipy.optimize 参数的错误
【发布时间】：2017-09-21 11:02:07
【问题描述】：

我将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 值正确包含拟合参数。如何计算与这些参数相关的误差？

【问题讨论】：

标签： numpy scipy curve-fitting

【解决方案1】：

这真的取决于你所说的“错误”是什么意思。您的问题没有一般性的答案，因为这取决于您适合什么以及您所做的假设。

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

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

【讨论】：

【解决方案2】：

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

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

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

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

【讨论】：

【解决方案3】：

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

这个方法的documentation 说

当(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol 时迭代停止。

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

当然，您正在应用像这样的最小化例程这一事实假设您的函数表现良好，从某种意义上说，它相当平滑并且找到的最优值很好地逼近接近最优值 em> 由参数 xⁱ 的二次函数：

其中Δxⁱ是参数xⁱ求得值与其最优值之差, H_ij 是 Hessian 矩阵。一点（令人惊讶的非平凡）线性代数可以让您获得一个非常标准的结果，以估计任何数量 X 的不确定性，这是您的参数 xⁱ 的函数：

这让我们可以写

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

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

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，您显然需要在此处进行更改。

【讨论】：