orthogonal distance from fit with scipy.odr

Question

I perform an ODR fit for given data points with x and y uncertainty with scipy.odr. For the fit result I would like to see the orthogonal distances of each point from the fit. This can be easily computed for a linear fit function but it can be more complicated for an arbitrary function. Since scipy.odr is minimizing these distances it has to compute them, so maybe one can directly access them but I did not find it in the documentation. Maybe I missed it? Any ideas?

So for the attached example code I would like to get an array of 100 elements containing the orthogonal distances and then e.g. plot a histogram.

import numpy as np
from scipy.odr import *

np.random.seed(1)

N = 100
x = np.linspace(0, 10, N)
y = 3*x - 1 + np.random.random(N)
sx = np.random.random(N)
sy = np.random.random(N)

def f(B, x):
    return B[0]*x + B[1]

linear = Model(f)
mydata = RealData(x, y, sx=sx, sy=sy)

myodr = ODR(mydata, linear, beta0=[1., 2.])
myoutput = myodr.run()
myoutput.pprint()

user_12750353 · Answer

Having a quick look in the documentation I think the model estimates for the orthogonal distance are what it brings input errors in myoutput.delta and output errors in myoutput.eps, combine the two of them we can get an estimate of the orthogonal distance.

Using the calculated y error, for this linear model we can use some geometry to find the orthogonal distance that is dy * sqrt(10/81)

dy = (y - f(myoutput.beta, x))
orthogonal = np.sqrt((dy / 3)**2 + (dy / 9)**2)
plt.plot(dy, np.sqrt(myoutput.delta**2 + myoutput.eps**2), '.', label='Estimated orthogonal distance')
plt.plot(dy, orthogonal, '.', label='Exact orthogonal distance')
plt.xlabel('Vertical distance')
plt.ylabel('Orthogonal distance')
plt.legend()

It is reasonable that the estimates are higher than the actual value, since it may be not estimating properly the closest point in the fitted curve.

orthogonal distance from fit with scipy.odr

1 Answers