1
2
3
4
5
6
7
8
9
10
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
36
37
38
39
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
47
48
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
128
129
130
132 r = 0.0
133 for i in range(len(params)):
134 r = r + params[i]*N.power(t, i)
135 return r
136
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
156
157 if __name__ == '__main__':
158
159 from Scientific.N import exp
160
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
172 return 1e13*exp(-param[0]/t)
173
174 print leastSquaresFit(f2, (3000.,), data_quantum)
175