python - How to compute standard error from ODR results?

Question

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]

Answer 1

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) = ?σ?2V β(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.