如何根据 ODR 结果计算标准误差？

Question

我使用 scipy.odr 是为了在这个问题之后对 x 和 y 的不确定性进行拟合 Correct fitting with scipy curve_fit including errors in x?

拟合后我想计算参数的不确定性。因此，我查看了协方差矩阵对角线元素的平方根。我得到：

>>> print(np.sqrt(np.diag(output.cov_beta)))
[ 0.17516591  0.33020487  0.27856021]

但是在 Output 中还有 output.sd_beta ，根据 doc on odr

Standard errors of the estimated parameters, of shape (p,).

但是，它并没有给我相同的结果：

>>> print(output.sd_beta)
[ 0.19705029  0.37145907  0.31336217]

编辑

这是笔记本上的示例：https://nbviewer.jupyter.org/github/gvallverdu/cookbook/blob/master/fit_odr.ipynb

最小二乘法

stop reason: ['Sum of squares convergence']
        params: [ -1.94792946  11.03369235  -5.43265555]
          info: 1
       sd_beta: [ 0.26176284  0.49877962  0.35510071]
sqrt(diag(cov): [ 0.25066236  0.47762805  0.34004208]

有 ODR

stop reason: ['Sum of squares convergence']
        params: [-1.93538595  6.141885   -3.80784384]
          info: 1
       sd_beta: [ 0.6941821   0.88909997  0.17292514]
sqrt(diag(cov): [ 0.01093697  0.01400794  0.00272447]

Answer 1

差异的原因是 sd_beta 被残差方差缩放，而 cov_beta 不是。

scipy.odr 是 ODRPACK FORTRAN 库的接口，它被薄包裹在 __odrpack.c 中。 sd_beta 和 cov_beta 通过索引到 FORTRAN 例程内部使用的 work 向量来恢复。它们在 work 中的第一个元素的索引是名为 sd 和 vcv 的变量（参见 here）。

来自 ODRPACK documentation（第 85 页）：

WORK(SDI) is the first element of a p × 1 array SD containing the standard deviations ̂σβK of the function parameters β, i.e., the square roots of the diagonal entries of the covariance matrix, where
WORK(SDI-1+K) = SD(K) = ̂V 1/2 β (K, K) = ̂σβK
for K = 1,... ,p.

WORK(VCVI) is the first element of a p × p array VCV containing the values of the covariance matrix of the parameters β prior to scaling by the residual variance, where
WORK(VCVI-1+I+(J-1)*(NP)) = VCV(I,J) = ̂σ⁻²V β(I, J)
for I = 1,... ,p and J = 1,... ,p.

换句话说，np.sqrt(np.diag(output.cov_beta * output.res_var)) 与 output.sd_beta 的结果相同。

我打开了错误报告 here。

如何根据 ODR 结果计算标准误差？

How to compute standard error from ODR results?

python

numpy

scipy

regression

curve-fitting