1 # Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
3 # This program is free software: you can redistribute it and/or modify
4 # it under the terms of the GNU General Public License as published by
5 # the Free Software Foundation, either version 3 of the License, or
6 # (at your option) any later version.
8 # This program is distributed in the hope that it will be useful,
9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 # GNU General Public License for more details.
13 # You should have received a copy of the GNU General Public License
14 # along with this program. If not, see <http://www.gnu.org/licenses/>.
16 # The author may be contacted at <wking@drexel.edu> on the Internet, or
17 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
18 # Philadelphia PA 19104, USA.
20 """Provide :class:`ModelFitter` to make arbitrary model fitting easy.
23 from copy import deepcopy
25 from numpy import arange, array, ndarray, exp, pi, sqrt
26 from scipy.optimize import leastsq
29 class ModelFitter (object):
30 """A convenient wrapper around :func:`scipy.optimize.leastsq`.
32 For minimizing fit error: returns the normalized rms difference
33 between a dataset and a function `y(X)`. In order to use this
34 class, you must sublcass it and overwrite `ModelFitter.model()`.
35 It is recommended that you also overwrite
36 `ModelFitter.guess_initial_params()`, and then not pass
37 initializing params to `ModelFitter.model()`. This helps ensure
38 good automation for unsupervised testing. Developing good
39 heuristics for `ModelFitter.guess_initial_params()` is usually the
40 hardest part writing a theoretical histogram test.
42 For bonus points, you can override `ModelFitter.print_model()` and
43 `ModelFitter.print_params()`, to give your users an easy to
44 understand idea of what's going on. It's better to be verbose and
45 clear. Remember, this is the testing code. It needs to be solid.
46 Save your wonderful, sneaky hacks for the main code ;).
48 TODO: look into :mod:`scipy.odr` as an alternative fitting
49 algorithm (minimizing the sum of squares of orthogonal distances,
50 vs. minimizing y distances).
55 Deflection data to be analyzed for the contact position.
57 Store any extra information useful inside your overridden
60 Rescale parameters so the guess for each is 1.0. Also rescale
61 the data so data.std() == 1.0.
65 >>> from pprint import pprint
67 >>> class LinearModel (ModelFitter):
68 ... '''Simple linear model.
70 ... Levenberg-Marquardt is not how you want to solve this problem
71 ... for real systems, but it's a simple test case.
73 ... def model(self, params):
74 ... '''A linear model.
78 ... .. math:: y = p_0 x + p_1
80 ... p = params # convenient alias
81 ... self._model_data[:] = p[0]*arange(len(self._data)) + p[1]
82 ... return self._model_data
83 ... def guess_initial_params(self):
84 ... return [float(self._data[-1] - self._data[0])/len(self._data),
86 ... def guess_scale(self, params):
88 ... if params[1] == 0:
91 ... offset_scale = 0.1*self._data.std()/abs(params[1])
92 ... if offset_scale == 0: # data is completely flat
94 ... return [slope_scale, offset_scale]
95 ... def model_string(self):
96 ... return 'y = A x + B'
97 ... def param_string(self, params):
99 ... for name,param in zip(['A', 'B'], params):
100 ... pstrings.append('%s=%.3f' % (name, param))
101 ... return ', '.join(pstrings)
102 >>> data = 20*numpy.sin(arange(1000)) + 7.*arange(1000) - 33.0
103 >>> m = LinearModel(data)
104 >>> slope,offset = m.fit()
105 >>> pprint(m.fit_info)
106 ... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE, +REPORT_UDIFF
107 {'active fitted parameters': array([ 6.999..., -32.889...]),
108 'active parameters': array([ 6.999..., -32.889...]),
109 'convergence flag': ...,
110 'covariance matrix': array([[ 1.199...e-08, -5.993...e-06],
111 [ -5.993...e-06, 3.994...e-03]]),
112 'data scale factor': 1.0,
113 'fitted parameters': array([ 6.999..., -32.889...]),
114 'info': {'fjac': array([[...]]),
115 'fvec': array([...]),
116 'ipvt': array([1, 2]),
118 'qtf': array([...])},
119 'initial parameters': [6.992..., -33.0],
120 'message': '...relative error between two consecutive iterates is at most 0.000...',
122 'scale': [0.100..., 6.123...]}
124 We round the outputs to protect the doctest against differences in
125 machine rounding during computation. We expect the values to be close
126 to the input settings (slope 7, offset -33).
128 >>> print '%.3f' % slope
130 >>> print '%.3f' % offset
133 >>> print m.model_string()
135 >>> print m.param_string([slope,offset])
137 >>> print ModelFitter.param_string(m, [slope,offset]) # doctest: +ELLIPSIS
138 param_0=6.9..., param_1=-32.8...
140 The offset is a bit off because, the range is not a multiple of
143 We could also use rescaled parameters:
145 >>> m = LinearModel(data, rescale=True)
146 >>> slope,offset = m.fit()
147 >>> print '%.3f' % slope
149 >>> print '%.3f' % offset
152 def __init__(self, *args, **kwargs):
153 self.set_data(*args, **kwargs)
155 def set_data(self, data, info=None, rescale=False):
157 self._model_data = ndarray(shape=data.shape, dtype=data.dtype)
159 self._rescale = rescale
161 for x in [data.std(), data.max()-data.min(), abs(data.max()), 1.0]:
163 self._data_scale_factor = x
166 self._data_scale_factor = 1.0
168 def model(self, params):
169 p = params # convenient alias
170 self._model_data[:] = arange(len(self._data))
171 raise NotImplementedError
173 def model_string(self):
174 """OVERRIDE: Give the user some idea of what model we're using.
176 return 'model string not defined'
178 def param_string(self, params):
179 """OVERRIDE: Give the user nicely named format for a parameter vector.
182 for i,p in enumerate(params):
183 pstrings.append('param_%d=%g' % (i,p))
184 return ', '.join(pstrings)
186 def guess_initial_params(self):
187 """OVERRIDE: Provide some intelligent heuristic to pick a
188 starting point for minimizing your model"""
189 raise NotImplementedError
191 def guess_scale(self, params):
192 """Guess the problem length scale in each parameter dimension.
196 From the :func:`scipy.optimize.leastsq` documentation, `diag`
197 (which we refer to as `scale`, sets `factor * || diag * x||`
198 as the initial step. If `x == 0`, then `factor` is used
199 instead (from `lmdif.f`_)::
202 wa3(j) = diag(j)*x(j)
206 if (delta .eq. zero) delta = factor
208 For most situations then, you don't need to do anything fancy.
209 The default scaling (if you don't set a scale) is::
211 c on the first iteration and if mode is 1, scale according
212 c to the norms of the columns of the initial jacobian.
214 (i.e. `diag(j) = acnorm(j)`, where `acnorm(j) is the norm of the `j`th column
215 of the initial Jacobian).
217 .. _lmdif.f: http://www.netlib.org/minpack/lmdif.f
221 def residual(self, params):
222 if self._rescale == True:
223 params = [p*s for p,s in zip(params, self._param_scale_factors)]
224 residual = self._data - self.model(params)
225 if False: # fit debugging
226 if not hasattr(self, '_i_'):
228 self._data.tofile('data.%d' % self._i_, sep='\n')
229 self._model_data.tofile('model.%d' % self._i_, sep='\n')
231 if self._rescale == True:
232 residual /= self._data_scale_factor
235 def fit(self, initial_params=None, scale=None, **kwargs):
239 initial_params : iterable or None
240 Initial parameter values for residual minimization. If
241 `None`, they are estimated from the data using
242 :meth:`guess_initial_params`.
243 scale : iterable or None
244 Parameter length scales for residual minimization. If
245 `None`, they are estimated from the data using
248 Any additional arguments are passed through to `leastsq`.
250 if initial_params == None:
251 initial_params = self.guess_initial_params()
253 scale = self.guess_scale(initial_params)
255 assert min(scale) > 0, scale
256 if self._rescale == True:
257 self._param_scale_factors = initial_params
258 for i,s in enumerate(self._param_scale_factors):
260 self._param_scale_factors[i] = 1.0
261 active_params = [p/s for p,s in zip(initial_params,
262 self._param_scale_factors)]
264 active_params = initial_params
265 params,cov,info,mesg,ier = leastsq(
266 func=self.residual, x0=active_params, full_output=True,
267 diag=scale, **kwargs)
268 if len(initial_params) == 1:
270 if self._rescale == True:
271 active_params = params
272 params = [p*s for p,s in zip(params, self._param_scale_factors)]
274 active_params = params
276 'rescaled': self._rescale,
277 'initial parameters': initial_params,
278 'active parameters': active_params,
280 'data scale factor': self._data_scale_factor,
281 'active fitted parameters': active_params,
282 'fitted parameters': params,
283 'covariance matrix': cov,
286 'convergence flag': ier,
291 class HistogramModelFitter (ModelFitter):
292 """Adapt ModelFitter for easy Histogram() comparison.
294 TODO: look into :func:`scipy.optimize.fmin_powell` as way to
295 minimize a single residul number. Useful for fitting `Histograms`
296 with their assorted residal methods.
300 >>> from random import gauss
301 >>> from numpy import array
302 >>> from .histogram import Histogram
304 >>> class GaussianModel (HistogramModelFitter):
305 ... '''Simple gaussian model.
307 ... def model(self, params):
308 ... '''A gaussian model.
312 ... .. math:: y \\propto e^{-\\frac{(x-p_0)^2}{2 p_1^2}}
314 ... p = params # convenient alias
315 ... p[1] = abs(p[1]) # cannot have negative std. dev.
316 ... self._model_data.counts = (
317 ... self.info['binwidth']*self.info['N']/(p[1]*sqrt(2*pi)) *
318 ... exp(-((self._model_data.bin_centers - p[0])/p[1])**2 / 2))
319 ... return self._model_data
320 ... def guess_initial_params(self):
321 ... return [self._data.mean, self._data.std_dev]
322 ... def guess_scale(self, params):
323 ... return [self._data.std_dev, self._data.std_dev]
324 ... def model_string(self):
325 ... return 'p(x) = A exp((x-mu)^2 / (2 sigma^2))'
326 ... def param_string(self, params):
328 ... for name,param in zip(['mu', 'sigma'], params):
329 ... pstrings.append('%s=%.3f' % (name, param))
330 ... return ', '.join(pstrings)
334 >>> data = array([gauss(mu,sigma) for i in range(int(1e4))])
336 >>> h.from_data(data, h.calculate_bin_edges(data, 1))
337 >>> m = GaussianModel(h)
339 >>> print m.model_string()
340 p(x) = A exp((x-mu)^2 / (2 sigma^2))
341 >>> print m.param_string(params) # doctest: +ELLIPSIS
342 mu=3.1..., sigma=1.4...
344 def set_data(self, data, info=None, rescale=False):
345 self._data = data # Histogram() instance
346 self._model_data = deepcopy(data)
347 self._model_data.bin_edges = array(self._model_data.bin_edges)
348 self._model_data.counts = array(self._model_data.counts)
352 self._rescale = rescale
353 self._data_scale_factor = 1.0 # bo need to rescale normalized hists.
354 self.info['N'] = self._data.total
355 # assume constant binwidth
356 self.info['binwidth'] = (
357 self._data.bin_edges[1] - self._data.bin_edges[0])
358 self._model_data.bin_centers = (
359 self._model_data.bin_edges[:-1] + self.info['binwidth']/2)
361 def residual(self, params):
362 if self._rescale == True:
363 params = [p*s for p,s in zip(params, self._param_scale_factors)]
364 residual = self._data.counts - self.model(params).counts
365 if False: # fit debugging
366 if not hasattr(self, '_i_'):
368 self._data.to_stream(open('hist-data.%d' % self._i_, 'w'))
369 self._model_data.headings = [str(x) for x in params]
370 self._model_data.to_stream(open('hist-model.%d' % self._i_, 'w'))
372 if self._rescale == True:
373 residual /= self._data_scale_factor