如何从ODR结果计算标准误差? [英] How to compute standard error from ODR results?

问题描述

为了解决这个问题，我使用scipy.odr来拟合x和y上的不确定性Output中也有output.sd_beta. html#scipy.odr.Output"rel =" noreferrer>关于odr的文档

形状(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]

解决方案

之所以出现差异，是因为sd_beta是根据残差进行缩放的，而cov_beta不是.

scipy.odr是ODRPACK FORTRAN库的接口，该库薄包装在__odrpack.c中.通过索引FORTRAN例程内部使用的work向量，可以恢复sd_beta和cov_beta. 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 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.

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屋！