Inconsistency in scipy.odr output for different fit types?

Question

For odr.set_job(fit_type=2) (ordinary least squares) odr.set_job(fit_type=0) (explicit ODR) the uncertainties sd_beta and the covariance matrix cov_beta are significantly different for the two fit types for the same data - which has an error in the y-part, only.

Has anyone an idea where this originates from? Below is a a code example for reproducing the effect. On the side: For the hidden scaling factor between sd_beta and cov_beta I found some help here: How to compute standard error from ODR results?

import numpy as np
import scipy.odr as spodr

def linear_func(pars, x):
    """Define a (linear) function to fit the data with."""
    a, b = pars
    return b * x + a

x = np.asfarray( [1.01585124, 2.03965699, 2.89148684] )
y = np.asfarray( [1.16201599, 2.15275179, 2.87248892] )
dy = np.asfarray( [0.15520856, 0.33160654, 0.45461028] )

real_data = spodr.RealData(x, y, sy=dy)
linear_model = spodr.Model(linear_func)
odr = spodr.ODR(real_data, linear_model, beta0=[0., 1.])
# fit_type=2 means least square fit
odr.set_job(fit_type=2)
out = odr.run()
print('fit_type=2 -------------------------------------------')
out.pprint()
print('\n')

# fit_type=0 means explicit ODR
odr.set_job(fit_type=0)
out = odr.run()
print('fit_type=0 -------------------------------------------')
out.pprint()

giving out the following result:

fit_type=2 -------------------------------------------
Beta: [0.22131257 0.92950773]
Beta Std Error: [0.0444006 0.03004  ]
Beta Covariance: [[ 0.10678885 -0.06586652]
 [-0.06586652  0.04888189]]
Residual Variance: 0.01846085213372428
Inverse Condition #: 0.02685695006039764
Reason(s) for Halting:
  Sum of squares convergence


fit_type=0 -------------------------------------------
Beta: [0.24510628 0.91688192]
Beta Std Error: [0.07210707 0.03482789]
Beta Covariance: [[ 2.29687687 -1.03022278]
 [-1.03022278  0.53584146]]
Residual Variance: 0.0022636956774060857
Inverse Condition #: 0.017439127336256896
Reason(s) for Halting:
  Sum of squares convergence

score 1 · Answer 1 · answered Jul 10 '20 at 08:52

It turned out that it is not a buggy inconsistency but more of a pitfall. For the ordinary least square type (OLS) of course the x-errors are effectively zero, because they just do not appear in the mathematical model.

Now, thinking that omitting sx in the ODR case implies sx → 0 is not correct. Actually for sx=None, sx is defaulted to unity, see here. This of course then yields something different compared to the OLS case.

Here, setting explicitly sx=data.x*ResEps as shown below gives the same result.

At minimum, this context should be noted in the docs somewhere more explicitly and I will try to prepare an issue on that. Interestingly in the code, the simultaneous passing of sx and covx raises a warning (which makes sense), but if neither of them is given, the above case occurs.

# obtain the machine precision for floating point operations
MachineEps = np.finfo('float').eps
ResEps = np.sqrt(MachineEps)

real_data = spodr.RealData(x, y, sy=dy, sx=x*ResEps)
linear_model = spodr.Model(linear_func)
odr = spodr.ODR(real_data, linear_model, beta0=[0., 1.])

# fit_type=0 means explicit ODR
odr.set_job(fit_type=0)
out = odr.run()
print('fit_type=0 -------------------------------------------')
out.pprint()

yields

fit_type=0 -------------------------------------------
Beta: [0.22131258 0.92950773]
Beta Std Error: [0.04440071 0.03004007]
Beta Covariance: [[ 0.10678941 -0.06586688]
 [-0.06586688  0.04888212]]
Residual Variance: 0.01846085213371618
Inverse Condition #: 0.02685687854396314
Reason(s) for Halting:
  Sum of squares convergence

Inconsistency in scipy.odr output for different fit types?

1 Answers1

Linked