Linear regression in 3 dimensions

[Robert Kern]

wirecom at wirelessmeasurement.com wrote:
> Hi all,
> 
> I am seeking a module that will do the equivalent of linear regression in
> 3D to yield a best fit a plane through a set of points (X1, Y1, Z1), (X1,
> Y1, Z1),... (Xn, Yn, Zn).
> 
> The resulting equation to be of the form:
> 
>  Z = aX + bY + c
> 
> The function I need would take the set of points and return a,c & c Any
> pointers to existing code / modules would be very helpful.

Well, that's a very unspecified problem. You haven't defined "best."

But if we make the assumption that you want to minimize the squared error in Z, 
that is minimize

   Sum((Z[i] - (a*X[i] + b*Y[i] + c)) ** 2)

then this is a standard linear algebra problem.

In [1]: import numpy as np

In [2]: a = 1.0

In [3]: b = 2.0

In [4]: c = 3.0

In [5]: rs = np.random.RandomState(1234567890)  # Specify a seed for repeatability

In [6]: x = rs.uniform(size=100)

In [7]: y = rs.uniform(size=100)

In [8]: e = rs.standard_normal(size=100)

In [9]: z = a*x + b*y + c + e

In [10]: A = np.column_stack([x, y, np.ones_like(x)])

In [11]: np.linalg.lstsq?
Type:           function
Base Class:     <type 'function'>
String Form:    <function lstsq at 0x6df070>
Namespace:      Interactive
File: 
/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/numpy-1.0b2.dev3002-py2.4-macosx-10.4-ppc.egg/numpy/linalg/linalg.py
Definition:     np.linalg.lstsq(a, b, rcond=1e-10)
Docstring:
     returns x,resids,rank,s
     where x minimizes 2-norm(|b - Ax|)
           resids is the sum square residuals
           rank is the rank of A
           s is the rank of the singular values of A in descending order

     If b is a matrix then x is also a matrix with corresponding columns.
     If the rank of A is less than the number of columns of A or greater than
     the number of rows, then residuals will be returned as an empty array
     otherwise resids = sum((b-dot(A,x)**2).
     Singular values less than s[0]*rcond are treated as zero.


In [12]: abc, residuals, rank, s = np.linalg.lstsq(A, z)

In [13]: abc
Out[13]: array([ 0.93104714,  1.96780364,  3.15185125])

-- 
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