Jak obliczyć błąd standardowy z wyników ODR?

Question

używam scipy.odr w celu dokonania dopasowanie z niepewnością na obu xiy po tym pytaniu Correct fitting with scipy curve_fit including errors in x?

Po pasowaniu chciałbym obliczyć niepewność na temat parametrów. Tak więc patrzę na pierwiastek kwadratowy z przekątnych elementów macierzy kowariancji. I otrzymujemy:

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

Ale w Output istnieje również output.sd_beta który jest, według doc on odr

błędy standardowe szacowanych parametrów, kształtu (P,).

Ale to nie daje mi takie same wyniki:

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

EDIT

To jest przykład na notebooku: https://nbviewer.jupyter.org/github/gvallverdu/cookbook/blob/master/fit_odr.ipynb

z najmniejszym placu

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] 

Z 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

Powodem rozbieżności jest to, że sd_beta jest skalowany przez resztkowej wariancji, natomiast cov_beta nie jest.

scipy.odr to interfejs biblioteki ODRPACK FORTRAN, który jest cienko zawinięty w __odrpack.c. sd_beta i cov_beta są odzyskiwane przez indeksowanie do wektora work, który jest używany wewnętrznie przez procedurę FORTRAN. Indeksy ich pierwszych elementów w work są zmiennymi o nazwach sd i vcv (patrz here).

Z ODRPACK documentation (str.85):

WORK(SDI) jest pierwszym elementem p × 1 tablicy SD zawierającego odchylenia standardowe ̂σβK funkcji Parametry β, tj kwadratowych korzenie przekątne wejścia macierzy kowariancji, gdzie

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

dla K = 1,... ,p.

WORK(VCVI) jest pierwszym elementem p × p tablicy VCV zawierającego wartości macierzy kowariancji parametrów βprzed skalowania od pozostałego wariancji gdzie

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

dla I = 1,... ,p i J = 1,... ,p.

Innymi słowy, np.sqrt(np.diag(output.cov_beta * output.res_var)) daje taki sam wynik, jak output.sd_beta.

Otworzyłem raport o błędzie here.