Source Code for Module gmisclib.s_lin_fit

 1  """Fit a plane to some data.""" 
 2   
 3   
 4  import g_lin_fit 
 5  import Num 
 6   
 7   
 8 -class Error(ValueError): 
 9 -        def __init__(self, s): 
10                  ValueError.__init__(self, s) 
11   
12   
13 -class NoDataError(ValueError): 
14 -        def __init__(self, s): 
15                  ValueError.__init__(self, s) 
16   
17 -def plane(data, wt=None): 
18          """Fit a plane to some data. 
19          where the fitting function is 
20          f = c[0]*1 + c[1]*dep[0] + c[2]*dep[1] + ... 
21          where c is the array of coefficients. 
22          The length of c is equal to the number of dependent 
23          variables in each datum. 
24          @param data: [ (independent, dependent, dependent...), ...] 
25          @return: (c, opt, resid, rank, sv), where 
26                  bestfit is the best fit to the data 
27                  (i.e., the values of f), 
28                  resid is the residual (a single float number), 
29                  rank the rank of the fit (int), 
30                  and sv is the array of the singular values. 
31          """ 
32          idata = [ ( [ t[0] for t in data ], None) ] 
33          n = len(data) 
34          if n <= 0: 
35                  raise NoDataError('No Data') 
36   
37          dims = len(data[0]) 
38          for t in data: 
39                  if len(t) != dims: 
40                          raise ValueError('Data length mismatch') 
41          ffi = Num.zeros((n, dims), Num.Float) 
42          ffi[:, 0] = 1 
43          for i in range(1, dims): 
44                  for j in range(n): 
45                          ffi[j, i] = data[j][i] 
46   
47          ff = [lambda info, ffin, nn, ii=i: ffin[:,ii] for i in range(dims) ] 
48   
49          opt, resid, rank, sv, bestfit = g_lin_fit.fit(idata, ff, ffi, wt) 
50          return (opt, bestfit[0], resid, rank, sv) 
51   
52   
53   
54 -def test(): 
55          data = [ (1, 0), (2, 2), (3, 4), (4, 6) ] 
56          opt, bestfit, resid, rank, sv = plane(data) 
57          assert rank == 2 
58          assert abs(resid) < 1e-6 
59          for i in range(len(data)): 
60                  assert abs(bestfit[i] - data[i][0]) < 1e-6 
61          assert abs(opt[0] - 1) < 1e-6 and abs(opt[1]-0.5) < 1e-6 
62   
63  if __name__ == '__main__': 
64          test() 
65