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

[Robert Kern]

Angus McMorland wrote:
> 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)?

It looks like you are doing the orthogonal distance version that I alluded to.
It depends on how you are measuring "error". If all of the variables ("input"
and "output") are all more or less the same kidn of "thing", and they together
form a Euclidean space with an isotropic metric (i.e., distance can be measured
simply by sqrt(sum(v*v)) for all vectors, with no weighting), then an
appropriate error function is the one you give: the shortest distance between
the point and the plane. Statistically speaking, you can frequently apply
weights and certain kinds of transformations to your variables to get this
condition. Assuming Gaussian errors (otherwise we shouldn't be using least
squares), you can derive the appropriate transformations from the covariance
matrices of the uncertainties.

In some models, however, the "input" variables are significantly different from
the "output" variables. Sometimes, they are simply incommensurate and you can't
really shove them together to form a real vector space with a metric. Other
times, you *could*, but the uncertainties on the "output" variables are so much
larger than those on the "input" variables that the only major contribution to
the residual is the difference in the "output" variables. In that case, the
classic "ordinary least squares" approach outlined in this thread using
linalg.lstsq() is appropriate.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco