[SciPy-User] Orthogonal distance regression in 3D

Robert Kern robert.kern@gmail....
Mon Mar 5 04:58:46 CST 2012


On Mon, Mar 5, 2012 at 10:26, Владимир <draft2008@bk.ru> wrote:
> 02 марта 2012, 15:49 от Robert Kern <robert.kern@gmail.com>:
>> On Fri, Mar 2, 2012 at 06:02, Владимир <draft2008@bk.ru> wrote:
>> > Hello!
>> > I'm working with orthogonal distance regression (scipy.odr).
>> > I try to fit the curve to a point cloud (3d), but it doesn work properly, it
>> > returns wrong results
>> >
>> > For example I want to fit the simple curve y = a*x + b*z + c to some point
>> > cloud (y_data, x_data, z_data)
>> >
>> >
>> >     def func(p, input):
>> >
>> >     x,z = input
>> >
>> >     x = np.array(x)
>> >
>> >     z = np.array(z)
>> >
>> >     return (p[0]*x + p[1]*z + p[2])
>> >
>> >
>> >     initialGuess = [1,1,1]
>> >
>> >     myModel = Model(func)
>> >
>> >     myData = Data([x_data, z_daya], y_data)
>> >
>> >     myOdr = ODR(myData, myModel, beta0 = initialGuess)
>> >
>> >     myOdr.set_job(fit_type=0)
>> >
>> >     out = myOdr.run()
>> >
>> >     print out.beta
>> >
>> > It works perfectly in 2d dimension (2 axes), but in 3d dimension the results
>> > are not even close to real, moreover it is very sensitive to initial Guess,
>> > so it returns different result even if i change InitiaGuess from [1,1,1]
>> > to [0.99,1,1]
>> >
>> > What do I do wrong?
>>
>> Can you provide a complete runnable example including some data? Note
>> that if you do not specify any errors on your data, they are assumed
>> to correspond to a standard deviation of 1 for all dimensions. If that
>> is wildly different from the actual variance around the "true"
>> surface, then it might lead the optimizer astray.
>>
>> --
>> Robert Kern
>>
>
> I wonder why when I change the initial guess the results changes too. As it, the result depends on the initial guess directly. This is wrong.
>
> Here is an example (Sorry for the huge array of data, but its important to show what happens on it)
>
> import numpy as np
> from scipy.odr import *
> from math import *

[snip]

> def funcReturner(p, input):
>        input = np.array(input)
>        x = input[0]
>        z = input[1]
>        return 10**(p[0]*x + p[1]*z +p[2])

Ah. 10**(p[0]*x+p[1]*z+p[2]) is a *lot* different from the linear
problem you initially asked about. Setting the uncertainties
accurately on all axes of your data is essential. Do you really know
what they are? It's possible that you want to try fitting a plane to
np.log10(y_data) instead.

> myModel = Model(funcReturner)
> myData = Data([x_data,z_data], y_data)
> myOdr = ODR(myData, myModel, beta0=[0.04, -0.02,  1.75])
> myOdr.set_job(fit_type=0)
> out = myOdr.run()
> result = out.beta
>
> print "Optimal coefficients: ", result
>
> I tryed to specify sx, sy, we, wd, delta, everything: and I get the better results, but they are still not what I need. And they are still depends directly on initial guess as well.
> If I set initial guess to [1,1,1], it fails to find any close solution and returns totally wrong result with huge Residual Variance like 3.21014784829e+209

For such a nonlinear problem, finding reasonable initial guesses is
useful. There is also a maximum iteration limit defaulting to a fairly
low 50. Check out.stopreason to see if it actually converged or just
ran into the iteration limit. You can keep calling myOdr.restart()
until it converges. If I start with beta0=[1,1,1], it converges
somewhere between 300 and 400 iterations.

-- 
Robert Kern


More information about the SciPy-User mailing list