Linear regression in 3 dimensions

[Andrew McLean]

Bernhard,

Levenberg-Marquardt is a good solution when you want to solve a general 
non-linear least-squares problem. As Robert said, the OPs problem is 
linear and Robert's solution exploits that. Using LM here is unnecessary 
and I suspect a fair bit less efficient (i.e. slower).

- Andrew


bernhard.voigt at gmail.com wrote:
> Hi Robert,
> 
> I'm using the scipy package for such problems. In the submodule
> scipy.optimize there is an implmentation of a least-square fitting
> algorithm (Levenberg-Marquardt) called leastsq.
> 
> You have to define a function that computes the residuals between your
> model and the data points:
> 
> import scipy.optimize
> 
> def model(parameter, x, y):
>     a, b, c = parameter
>     return a*x + b*y + c
> 
> def residual(parameter, data, x, y):
>     res = []
>     for _x in x:
>         for _y in y:
>             res.append(data-model(parameter,x,y)
>     return res
> 
> params0 = [1., 1.,1.]
> result = scipy.optimize.leastsq(resdiual, params0, (data,x,y))
> fittedParams = result[0]
> 
> If you haven't used numeric, numpy or scipy before, you should take a
> look at an introduction. It uses some nice extended array objects,
> where you can use some neat index tricks and compute values of array
> items without looping through it.
> 
> Cheers! Bernhard
> 
> 
> 
> Robert Kern wrote:
>> 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
>