Linear least-squares estimator in Python

[Will Ware]

import Numeric, LinearAlgebra

class LLSExyz:
"""Linear least-square estimator mapping points (x,y) to result z
"""
def __init__(self, xylist, zlist):
assert len(xylist) == len(zlist)
n = len(xylist)
sx = sy = sz = 0.
sxx = syy = sxy = sxz = syz = 0.
for i in range(n):
x, y = xylist[i]
z = zlist[i]
sx += x
sy += y
sz += z
sxx += x * x
syy += y * y
sxy += x * y
sxz += x * z
syz += y * z
A = Numeric.array(((sxx, sxy, sx),
(sxy, syy, sy),
(sx, sy, n)), Numeric.Float)
Ainv = LinearAlgebra.inverse(A)
B = Numeric.array(((sxz,),
(syz,),
(sz,)), Numeric.Float)
abc = Numeric.matrixmultiply(Ainv, B)
self.a, self.b, self.c = map(lambda x: x[0], abc)
def abc(self):
return (self.a, self.b, self.c)
def __call__(self, x, y):
return self.a * x + self.b * y + self.c

#########

gpsData = (
(42.53071, -71.28185),
(42.52875, -71.27457),
(42.53498, -71.26556),
(42.53200, -71.25793),
(42.53752, -71.25861),
(42.54243, -71.25229),
(42.53411, -71.23232)
)

pixelData = (
(46, 235),
(108, 246),
(215, 147),
(304, 192),
(294, 107),
(369, 28),
(601, 166)
)

f = LLSExyz(gpsData, map(lambda x: x[0], pixelData))
g = LLSExyz(gpsData, map(lambda x: x[1], pixelData))

for i in range(len(gpsData)):
x, y = gpsData[i]
lat, long = pixelData[i]
print lat, f(x, y), long, g(x, y)