如何从 ODR 结果计算标准误差? [英] How to compute standard error from ODR results?
问题描述
我使用 scipy.odr
以适应 x 和 y 上的不确定性,并遵循这个问题 使用 scipy curve_fit 进行正确拟合,包括 x 中的误差?
拟合后,我想计算参数的不确定性.因此,我查看协方差矩阵的对角元素的平方根.我明白了:
<预><代码>>>>打印(np.sqrt(np.diag(output.cov_beta)))[ 0.17516591 0.33020487 0.27856021]但是在 Output
中还有 output.sd_beta
,根据 odr 上的文档
估计参数的标准误差,形状为 (p,).
但是,它并没有给我相同的结果:
<预><代码>>>>打印(输出.sd_beta)[ 0.19705029 0.37145907 0.31336217]编辑
这是笔记本上的示例:https://nbviewer.jupyter.org/github/gvallverdu/cookbook/blob/master/fit_odr.ipynb
最小二乘
停止原因:['平方和收敛']参数:[-1.94792946 11.03369235 -5.43265555]信息:1sd_beta:[0.26176284 0.49877962 0.35510071]sqrt(diag(cov): [ 0.25066236 0.47762805 0.34004208]
有 ODR
停止原因:['平方和收敛']参数:[-1.93538595 6.141885 -3.80784384]信息:1sd_beta:[0.6941821 0.88909997 0.17292514]sqrt(diag(cov): [ 0.01093697 0.01400794 0.00272447]
产生差异的原因是 sd_beta
是由残差方差缩放的,而 cov_beta
不是'
scipy.odr
是 ODRPACK FORTRAN 库的接口,它被薄薄地包裹在 __odrpack.c
中.sd_beta
和 cov_beta
通过索引到 FORTRAN 例程内部使用的 work
向量来恢复.work
中第一个元素的索引是名为 sd
和 vcv
的变量(参见 这里).
来自 ODRPACK 文档(第 85 页):
<块引用>WORK(SDI)
是 p × 1
数组 SD
的第一个元素,包含函数参数β
的标准差̂σβK
,即协方差矩阵对角线项的平方根,其中
WORK(SDI-1+K) = SD(K) = ̂V 1/2 β (K, K) = ̂σβK
对于 K = 1,... ,p
.
WORK(VCVI)
是 p × p
数组 VCV
的第一个元素,包含参数β
之前的协方差矩阵的值通过残差进行缩放,其中
WORK(VCVI-1+I+(J-1)*(NP)) = VCV(I,J) = ̂σ⁻²V β(I, J)
对于 I = 1,... ,p
和 J = 1,... ,p
.
换句话说,np.sqrt(np.diag(output.cov_beta * output.res_var))
会给你与 output.sd_beta
相同的结果.
我在此处打开了错误报告.
I use scipy.odr
in order to make a fit with uncertainties on both x and y following this question Correct fitting with scipy curve_fit including errors in x?
After the fit I would like to compute the uncertainties on the parameters. Thus I look at the square root of the diagonal elements of the covariance matrix. I get :
>>> print(np.sqrt(np.diag(output.cov_beta)))
[ 0.17516591 0.33020487 0.27856021]
But in the Output
there is also output.sd_beta
which is, according to the doc on odr
Standard errors of the estimated parameters, of shape (p,).
But, it does not give me the same results :
>>> print(output.sd_beta)
[ 0.19705029 0.37145907 0.31336217]
EDIT
This is an example on a notebook : https://nbviewer.jupyter.org/github/gvallverdu/cookbook/blob/master/fit_odr.ipynb
With least square
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]
With 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]
The reason for the discrepancy is that sd_beta
is scaled by the residual variance, whereas cov_beta
isn't.
scipy.odr
is an interface for the ODRPACK FORTRAN library, which is thinly wrapped in __odrpack.c
. sd_beta
and cov_beta
are recovered by indexing into the work
vector that's used internally by the FORTRAN routine. The indices of their first elements in work
are variables named sd
and vcv
(see here).
From the ODRPACK documentation (p.85):
WORK(SDI)
is the first element of ap × 1
arraySD
containing the standard deviationŝσβK
of the function parametersβ
, i.e., the square roots of the diagonal entries of the covariance matrix, whereWORK(SDI-1+K) = SD(K) = ̂V 1/2 β (K, K) = ̂σβK
for
K = 1,... ,p
.
WORK(VCVI)
is the first element of ap × p
arrayVCV
containing the values of the covariance matrix of the parametersβ
prior to scaling by the residual variance, whereWORK(VCVI-1+I+(J-1)*(NP)) = VCV(I,J) = ̂σ⁻²V β(I, J)
for
I = 1,... ,p
andJ = 1,... ,p
.
In other words, np.sqrt(np.diag(output.cov_beta * output.res_var))
will give you the same result as output.sd_beta
.
I've opened a bug report here.
这篇关于如何从 ODR 结果计算标准误差?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!