[SciPy-user] using linalg.lstsq (and now optimize.leastsq)

[Angus McMorland]

Hi all,

Bill Baxter wrote:
> Try this:
> 
> #                x   y
> A = array([[ 0., 0.],
>                 [ 0., 1.],
>                 [ 1., 0.],
>                 [ 1., 1.]])
> A = hstack( (A, ones((A.shape[0],1)) ))
> z = array([.5,.6,.4,.55])
> c = scipy.linalg.lstsq(A,z)[0]
> 

<snip>

>     Mitchell Timin wrote:
> 
>     > I accept that, but I don't see how to put my problem into that
>     form.  Suppose we are given these 4 points only:
>     >    x     y      z
>     >    0     0      .5
>     >    0     1      .6
>     >    1     0      .4
>     >    1     1      .55
>     >

Just for fun, I tried to do this using optimize.leastsq as follows (a
method I find conceptually simpler to follow):

from numpy import *
from scipy.optimize import leastsq

def resids( p, pts ):
    '''returns the residual (distances to the plane for the
    points in pts) for a plane given by ax + by + cz + d = 0'''

    (a, b, c, d) = p
    return (a * pts[:,0] + b* pts[:,1] + c * pts[:,2] + d) \
           / sqrt(a**2 + b**2 + c**2)

def fit_plane( pts, p0=None ):
    if p0 == None:
        p0 = (1,1,1,1)
    plsq = leastsq( resids, p0, args=(pts,) )
    return plsq

then:

pts = array([[0,0,0.5],[0,1,0.6],[1,0,0.4],[1,1,0.55]])
plfit = fit_plane( pts )[0]
print -plfit / plfit[2]   # convert to the other form of plane eqn
-> [-0.07504592  0.12507655 -1.         0.48748469]

print c # from Bill's code
-> [-0.075   0.125   0.4875]

Is the difference between the two answers indicative of anything wrong
(or does it suggest the better applicability of one method over the other)?

Cheers,

Angus.
-- 
Angus McMorland
email a.mcmorland at auckland.ac.nz
mobile +64-21-155-4906

PhD Student, Neurophysiology / Multiphoton & Confocal Imaging
Physiology, University of Auckland
phone +64-9-3737-599 x89707

Armourer, Auckland University Fencing
Secretary, Fencing North Inc.