CODEWITHSUNDEEP
pythonscipyregressioncurve-fittingorthogonal
Chryssi Koukouraki
Chryssi Koukouraki

Reputation: 31

perpendicular residuals with scipy odr

I want to fit data to a power-law using scipy odr, which aims to minimize the perpendicular squared distances. I want to find the perpendicular distance between each point and the fit.

The ".delta" returns the |y_model - y_data| residual (vertical distances) for each point and I can't seem to find any built-in method that returns the perpendicular distances. I am sure that they are calculated internally in the algorithm, because they are used to calculate the parameters of the fit, but I am not able to find them anywhere.

If there is indeed no way to get them with one command, how can I do this manually? I can do it if the fitted curve is a straight line, but it is a power-law, and I don't know what to do.

Upvotes: 3

Views: 100

Answers (1)

jlandercy
jlandercy

Reputation: 11002

Lets say you have the following model:

import numpy as np
import matplotlib.pyplot as plt
from scipy import odr
from shapely.geometry import Point, LineString

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

np.random.seed(12345)
x = np.linspace(0, 1, 20)
y = model([0.1, 0.2, 1.], x)
sx = 0.1 * np.ones_like(x)
sy = 0.2 * np.ones_like(x)
x += sx * np.random.normal(size=x.size)
y += sy * np.random.normal(size=x.size)

regressor = odr.Model(model)
data = odr.RealData(x, y, sx=sx, sy=sy)
solver = odr.ODR(data, regressor, beta0=[1., 1., 1.])
solution = solver.run()

Inspecting the solution object reveals no residuals (distances).

Then we can use shapely to finish the job.

First we resample and estimate the model (with a bit of extrapolation):

xlin = np.linspace(-0.2, 1.1, 200)
ylin = model(solution.beta, xlin)

We construct a LineString from the interpolation:

fit = LineString([(x, y) for x, y in zip(xlin, ylin)])

Now we can easily compute distance from this curve:

distances = np.array([
    fit.distance(Point(x, y))
    for x, y in zip(x, y)
])

Or even better locate point projections on the curve:

projected = [
    fit.interpolate(fit.project(Point(x, y)))
    for x, y in zip(x, y)
]
projected = np.array([[p.x, p.y] for p in projected])

Graphically it renders as:

enter image description here

Upvotes: 1

Related Questions