14

The leastsq method in scipy lib fits a curve to some data. And this method implies that in this data Y values depends on some X argument. And calculates the minimal distance between curve and the data point in the Y axis (dy)

But what if I need to calculate minimal distance in both axes (dy and dx)

Is there some ways to implement this calculation?

Here is a sample of code when using one axis calculation:

import numpy as np
from scipy.optimize import leastsq

xData = [some data...]
yData = [some data...]

def mFunc(p, x, y):
    return y - (p[0]*x**p[1])  # is takes into account only y axis

plsq, pcov = leastsq(mFunc, [1,1], args=(xData,yData))
print plsq

I recently tryed scipy.odr library and it returns the proper results only for linear function. For other functions like y=a*x^b it returns wrong results. This is how I use it:

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

myModel = Model(f)
myData = Data(xData, yData)
myOdr = ODR(myData, myModel , beta0=[1,1])
myOdr.set_job(fit_type=0) #if set fit_type=2, returns the same as leastsq
out = myOdr.run()
out.pprint()

This returns wrong results, not desired, and in some input data not even close to real. May be, there is some special ways of using it, what do I do wrong?

Vladimir
  • 601
  • 2
  • 7
  • 13
  • 1
    Scipy does have a module for "Orthogonal Distance Regression" - is that what you need? http://docs.scipy.org/doc/scipy/reference/odr.html – Thomas K Feb 21 '12 at 12:35
  • Yes, it seems to resolve this problem, but when I try it, it returns the same result as leastsq method. I followed the examples, which are given in documentation, and it doesnt work as needed. Do you have some working examples? – Vladimir Feb 22 '12 at 04:41
  • When I tried it I found that the results were similar, but not identical - I assumed that just meant that the extra calculation didn't make much difference to the fit. – Thomas K Feb 22 '12 at 12:33
  • 3
    I got It! I've found the solution. The problem was in inappropriate initial guesses for odr solver(beta0 parameter). – Vladimir Feb 24 '12 at 05:29

4 Answers4

9

I've found the solution. Scipy Odrpack works noramally but it needs a good initial guess for correct results. So I divided the process into two steps.

First step: find the initial guess by using ordinaty least squares method.

Second step: substitude these initial guess in ODR as beta0 parameter.

And it works very well with an acceptable speed.

Thank you guys, your advice directed me to the right solution

Vladimir
  • 601
  • 2
  • 7
  • 13
  • 3
    I know this is an old post, but could you possibly post your code snippet here. I am trying to to do an implicit ODR but I am not sure how to set it up in scipy. – Barbarossa Dec 06 '13 at 19:12
  • 1
    @Barbarossa You might like [this code snippet](http://blog.rtwilson.com/orthogonal-distance-regression-in-python/). – gerrit Aug 09 '16 at 20:54
  • 2
    This is not really an answer. You should/could have shared the code as well. – Foad S. Farimani Feb 27 '20 at 21:13
8

scipy.odr implements the Orthogonal Distance Regression. See the instructions for basic use in the docstring and documentation.

Foad S. Farimani
  • 12,396
  • 15
  • 78
  • 193
Robert Kern
  • 13,118
  • 3
  • 35
  • 32
  • oh, yeah, I tryed it, and it works the same way as leastsq, returns the same results – Vladimir Feb 22 '12 at 04:03
  • 4
    The exact same results? I find that unlikely. Can you update your post to show a runnable example with `scipy.odr` that gives you wrong results? – Robert Kern Feb 22 '12 at 11:43
0

you can use the ONLS package in R.

ly g
  • 13
  • 3
  • 1
    Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers [in the help center](/help/how-to-answer). – Community Apr 20 '22 at 07:51
0

If/when you are able to invert the function described by p you may just include x-pinverted(y) in mFunc, I guess as sqrt(a^2+b^2), so (pseudo code)

return sqrt( (y - (p[0]*x**p[1]))^2 + (x - (pinverted(y))^2)

for example for 

y=kx+m   p=[m,k]    
pinv=[-m/k,1/k]

return sqrt( (y - (p[0]+x*p[1]))^2 + (x - (pinv[0]+y*pinv[1]))^2)

But what you ask for is in some cases problematic. For example, if a polynomial (or your x^j) curve has a minimum ym at y(m) and you have a point x,y lower than ym, what kind of value do you want to return? There's not always a solution.

Johan Lundberg
  • 26,184
  • 12
  • 71
  • 97
  • If the point x,y lower than ym then it should return the miminal distance to that ym. **sqrt((m-x)^2 + (ym-y)^2)/2**. Why this is a problem? – Vladimir Feb 21 '12 at 12:33
  • if you have a function y=f(x) that means that for any x there is a value for y. Bot there's not always a value of x for any input y. Not all functions are invertible. for example y=x^2 and any point with x=-2 – Johan Lundberg Feb 21 '12 at 12:39
  • oh, I now I see, yes you are right. But how some software build this total least square (Deming regression) fitting for quite any function on any input data. There must be some way to do it in python – Vladimir Feb 21 '12 at 12:47
  • I think Deming regression is linear so inversion is not complicated. So, it's similar to the linear example I gave. Your curve is not linear an it's not invertible. – Johan Lundberg Feb 21 '12 at 12:52