¿Cómo calcular el error estándar de los resultados de ODR?

Uso scipy.odr para hacer un ajuste con incertidumbres tanto en x como en y después de esta pregunta ¿ Corregir el ajuste con scipy curve_fit incluyendo errores en x?

Después del ajuste, me gustaría calcular las incertidumbres sobre los parámetros. Así miro la raíz cuadrada de los elementos diagonales de la matriz de covarianza. Yo obtengo :

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

Pero en la Output también hay output.sd_beta que es, de acuerdo con el documento en odr

Errores estándar de los parámetros estimados, de forma (p,).

Pero, no me da los mismos resultados:

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

EDITAR

Este es un ejemplo en un cuaderno: https://nbviewer.jupyter.org/github/gvallverdu/cookbook/blob/master/fit_odr.ipynb

Con menos cuadrado

 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] 

Con 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] 

La razón de la discrepancia es que sd_beta es escalada por la varianza residual, mientras que cov_beta no lo es.

scipy.odr es una interfaz para la biblioteca ODRPACK FORTRAN, que está finamente envuelta en __odrpack.c . sd_beta y cov_beta se recuperan al indexar el vector de work que es utilizado internamente por la rutina FORTRAN. Los índices de sus primeros elementos en el work son variables denominadas sd y vcv (ver aquí ).

De la documentación de ODRPACK (p.85):

WORK(SDI) es el primer elemento de una matriz p × 1 SD contiene las desviaciones estándar ̂σβK de los parámetros de función β , es decir, las raíces cuadradas de las entradas diagonales de la matriz de covarianza, donde

 WORK(SDI-1+K) = SD(K) = ̂V 1/2 β (K, K) = ̂σβK 

para K = 1,... ,p .

WORK(VCVI) es el primer elemento de una matriz p × p VCV contiene los valores de la matriz de covarianza de los parámetros β antes de escalar por la varianza residual , donde

 WORK(VCVI-1+I+(J-1)*(NP)) = VCV(I,J) = ̂σ⁻²V β(I, J) 

para I = 1,... ,p y J = 1,... ,p .

En otras palabras, np.sqrt(np.diag(output.cov_beta * output.res_var)) le dará el mismo resultado que output.sd_beta .

He abierto un informe de error aquí .