[SciPy-user] using linalg.lstsq (and now optimize.leastsq)
Angus McMorland
a.mcmorland at auckland.ac.nz
Thu May 11 01:02:11 EDT 2006
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.
More information about the SciPy-User
mailing list