[frames] | no frames]

Source Code for Module nictools.SP_LeastSquares

1  # This module contains functions to do general non-linear
2  # least squares fits.
3  #
4  # Written by Konrad Hinsen <hinsen@cnrs-orleans.fr>
5  # last revision: 2006-11-23
6  #
7  # Modified import statements to use just a portion of the package
8  # and drop support for non-numpy array packages. Also added an
10  # V. Laidler, 2007-02-05.
11
12  """
13  Non-linear least squares fitting
14
15  Usage example::
16
17      from Scientific.N import exp
18
19      def f(param, t):
20          return param[0]*exp(-param[1]/t)
21
22      data_quantum = [(100, 3.445e+6),(200, 2.744e+7),
23                      (300, 2.592e+8),(400, 1.600e+9)]
24      data_classical = [(100, 4.999e-8),(200, 5.307e+2),
25                        (300, 1.289e+6),(400, 6.559e+7)]
26
27      print leastSquaresFit(f, (1e13,4700), data_classical)
28
29      def f2(param, t):
30          return 1e13*exp(-param[0]/t)
31
32      print leastSquaresFit(f2, (3000.,), data_quantum)
33  """
34  #------------------------------------------------------
35  #Original import statements
36  #------------------------------------------------------
37  #from Scientific import N, LA
38  #from FirstDerivatives import DerivVar
39  #from Scientific import IterationCountExceededError
40  #------------------------------------------------------
41  from __future__ import division
42  import SP_numpy as N
43  import numpy.oldnumeric.linear_algebra as LA
44  from SP_FirstDerivatives import DerivVar
45 -class IterationCountExceededError(ValueError):
46 pass
47 48
49 -def _chiSquare(model, parameters, data):
50 n_param = len(parameters) 51 chi_sq = 0. 52 alpha = N.zeros((n_param, n_param)) 53 for point in data: 54 sigma = 1 55 if len(point) == 3: 56 sigma = point[2] 57 f = model(parameters, point[0]) 58 chi_sq = chi_sq + ((f-point[1])/sigma)**2 59 d = N.array(f[1])/sigma 60 alpha = alpha + d[:,N.NewAxis]*d 61 return chi_sq, alpha
62
63 -def leastSquaresFit(model, parameters, data, max_iterations=None, 64 stopping_limit = 0.005):
65 """General non-linear least-squares fit using the 66 X{Levenberg-Marquardt} algorithm and X{automatic differentiation}. 67 68 @param model: the function to be fitted. It will be called 69 with two parameters: the first is a tuple containing all fit 70 parameters, and the second is the first element of a data point (see 71 below). The return value must be a number. Since automatic 72 differentiation is used to obtain the derivatives with respect to the 73 parameters, the function may only use the mathematical functions known 74 to the module FirstDerivatives. 75 @type param: callable 76 77 @param parameters: a tuple of initial values for the 78 fit parameters 79 @type parameters: C{tuple} of numbers 80 81 @param data: a list of data points to which the model 82 is to be fitted. Each data point is a tuple of length two or 83 three. Its first element specifies the independent variables 84 of the model. It is passed to the model function as its first 85 parameter, but not used in any other way. The second element 86 of each data point tuple is the number that the return value 87 of the model function is supposed to match as well as possible. 88 The third element (which defaults to 1.) is the statistical 89 variance of the data point, i.e. the inverse of its statistical 90 weight in the fitting procedure. 91 @type data: C{list} 92 93 @returns: a list containing the optimal parameter values 94 and the chi-squared value describing the quality of the fit 95 @rtype: C{(list, float)} 96 """ 97 n_param = len(parameters) 98 p = () 99 i = 0 100 for param in parameters: 101 p = p + (DerivVar(param, i),) 102 i = i + 1 103 id = N.identity(n_param) 104 l = 0.001 105 chi_sq, alpha = _chiSquare(model, p, data) 106 niter = 0 107 itertrace=[p] 108 while 1: 109 delta = LA.solve_linear_equations(alpha+l*N.diagonal(alpha)*id, 110 -0.5*N.array(chi_sq[1])) 111 next_p = map(lambda a,b: a+b, p, delta) 112 next_chi_sq, next_alpha = _chiSquare(model, next_p, data) 113 if next_chi_sq > chi_sq: 114 l = 10.*l 115 else: 116 l = 0.1*l 117 if chi_sq[0] - next_chi_sq[0] < stopping_limit: break 118 p = next_p 119 chi_sq = next_chi_sq 120 alpha = next_alpha 121 niter = niter + 1 122 if max_iterations is not None and niter == max_iterations: 123 raise IterationCountExceededError 124 return [p[0] for p in next_p], next_chi_sq[0], itertrace
125 126 # 127 # The special case of n-th order polynomial fits 128 # was contributed by David Ascher. Note: this could also be 129 # done with linear least squares, e.g. from LinearAlgebra. 130 #
131 -def _polynomialModel(params, t):
132 r = 0.0 133 for i in range(len(params)): 134 r = r + params[i]*N.power(t, i) 135 return r
136
137 -def polynomialLeastSquaresFit(parameters, data):
138 """ 139 Least-squares fit to a polynomial whose order is defined by 140 the number of parameter values. 141 142 @note: This could also be done with a linear least squares fit 143 from L{LinearAlgebra} 144 145 @param parameters: a tuple of initial values for the polynomial 146 coefficients 147 @type parameters: C{tuple} 148 149 @param data: the data points, as for L{leastSquaresFit} 150 @type data: C{list} 151 """ 152 return leastSquaresFit(_polynomialModel, parameters, data)
153 154 155 # Test code 156 157 if __name__ == '__main__': 158 159 from Scientific.N import exp 160
161 - def f(param, t):
162 return param[0]*exp(-param[1]/t)
163 164 data_quantum = [(100, 3.445e+6),(200, 2.744e+7), 165 (300, 2.592e+8),(400, 1.600e+9)] 166 data_classical = [(100, 4.999e-8),(200, 5.307e+2), 167 (300, 1.289e+6),(400, 6.559e+7)] 168 169 print leastSquaresFit(f, (1e13,4700), data_classical) 170
171 - def f2(param, t):
172 return 1e13*exp(-param[0]/t)
173 174 print leastSquaresFit(f2, (3000.,), data_quantum) 175

 Generated by Epydoc 3.0.1 on Thu Jan 13 11:20:16 2011 http://epydoc.sourceforge.net