Linear regression in 3 dimensions
andrew-news at andros.org.uk
Fri Sep 8 21:17:00 CEST 2006
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).
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:
> return res
> params0 = [1., 1.,1.]
> result = scipy.optimize.leastsq(resdiual, params0, (data,x,y))
> fittedParams = result
> 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 : import numpy as np
>> In : a = 1.0
>> In : b = 2.0
>> In : c = 3.0
>> In : rs = np.random.RandomState(1234567890) # Specify a seed for repeatability
>> In : x = rs.uniform(size=100)
>> In : y = rs.uniform(size=100)
>> In : e = rs.standard_normal(size=100)
>> In : z = a*x + b*y + c + e
>> In : A = np.column_stack([x, y, np.ones_like(x)])
>> In : np.linalg.lstsq?
>> Type: function
>> Base Class: <type 'function'>
>> String Form: <function lstsq at 0x6df070>
>> Namespace: Interactive
>> Definition: np.linalg.lstsq(a, b, rcond=1e-10)
>> 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*rcond are treated as zero.
>> In : abc, residuals, rank, s = np.linalg.lstsq(A, z)
>> In : abc
>> Out: 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
More information about the Python-list