Source code for numkit.fitting

# numkit --- data fitting
# Copyright (c) 2010 Oliver Beckstein <>
# Released under the "Modified BSD Licence" (see COPYING).

import numpy
import logging
logger = logging.getLogger("numkit.fitting")

[docs]def Pearson_r(x,y): """Pearson's r (correlation coefficient). Pearson(x,y) --> correlation coefficient *x* and *y* are arrays of same length. Historical note -- Naive implementation of Pearson's r :: Ex = scipy.stats.mean(x) Ey = scipy.stats.mean(y) covxy = numpy.sum((x-Ex)*(y-Ey)) r = covxy/math.sqrt(numpy.sum((x-Ex)**2)*numpy.sum((y-Ey)**2)) """ return numpy.corrcoef(x,y)[1,0]
[docs]def linfit(x,y,dy=None): """Fit a straight line y = a + bx to the data in *x* and *y*. Errors on y should be provided in dy in order to assess the goodness of the fit and derive errors on the parameters. linfit(x,y[,dy]) --> result_dict Fit y = a + bx to the data in x and y by analytically minimizing chi^2. dy holds the standard deviations of the individual y_i. If dy is not given, they are assumed to be constant (note that in this case Q is set to 1 and it is meaningless and chi2 is normalised to unit standard deviation on all points!). Returns the parameters a and b, their uncertainties sigma_a and sigma_b, and their correlation coefficient r_ab; it also returns the chi-squared statistic and the goodness-of-fit probability Q (that the fit would have chi^2 this large or larger; Q < 10^-2 indicates that the model is bad --- Q is the probability that a value of chi-square as _poor_ as the calculated statistic chi2 should occur by chance.) :Returns: result_dict with components intercept, sigma_intercept a +/- sigma_a slope, sigma_slope b +/- sigma_b parameter_correlation correlation coefficient r_ab between a and b chi_square chi^2 test statistic Q_fit goodness-of-fit probability Based on 'Numerical Recipes in C', Ch 15.2. """ if dy is None: dy = [] import scipy.stats n = len(x) m = len(y) if n != m: raise ValueError("lengths of x and y must match: {0!s} != {1!s}".format(n, m)) try: have_dy = (len(dy) > 0) except TypeError: have_dy = False if not have_dy: dy = numpy.ones((n), numpy.float64) x = numpy.asarray(x) y = numpy.asarray(y) dy = numpy.asarray(dy) s2 = dy*dy S = numpy.add.reduce(1/s2) Sx = numpy.add.reduce(x/s2) Sy = numpy.add.reduce(y/s2) Sxx = numpy.add.reduce(x*x/s2) Sxy = numpy.add.reduce(x*y/s2) t = (x - Sx/S)/dy Stt = numpy.add.reduce(t*t) b = numpy.add.reduce(t*y/dy)/Stt a = (Sy - Sx*b)/S sa = numpy.sqrt((1 + (Sx*Sx)/(S*Stt))/S) sb = numpy.sqrt(1/Stt) covab = -Sx/(S*Stt) r = covab/(sa*sb) chi2 = numpy.add.reduce(((y-a-b*x)/dy)**2) if not have_dy: # estimate error if none were provided sigmadata = numpy.sqrt(chi2/(n-2)) sa *= sigmadata sb *= sigmadata Q = 1.0 else: Q = scipy.stats.distributions.chi2.sf(chi2, n-2) return {"intercept":a,"slope":b, "sigma_intercept":sa,"sigma_slope":sb, "parameter_correlation":r, "chi_square":chi2, "Q":Q}
[docs]class FitFunc(object): """Fit a function f to data (x,y) using the method of least squares. The function is fitted when the object is created, using :func:`scipy.optimize.leastsq`. One must derive from the base class :class:`FitFunc` and override the :meth:`FitFunc.f_factory` (including the definition of an appropriate local :func:`fitfunc` function) and :meth:`FitFunc.initial_values` appropriately. See the examples for a linear fit :class:`FitLin`, a 1-parameter exponential fit :class:`FitExp`, or a 3-parameter double exponential fit :class:`FitExp2`. The object provides two attributes :attr:`FitFunc.parameters` list of parameters of the fit :attr:`FitFunc.message` message from :func:`scipy.optimize.leastsq` After a successful fit, the fitted function can be applied to any data (a 1D-numpy array) with :meth:``. """ def __init__(self,x,y,parameters=None): import scipy.optimize _x = numpy.asarray(x) _y = numpy.asarray(y) p0 = self._get_initial_values(parameters) fitfunc = self.f_factory() def errfunc(p,x,y): return fitfunc(p,x) - y # residuals p,msg = scipy.optimize.leastsq(errfunc,p0[:],args=(_x,_y)) try: p[0] self.parameters = p except (TypeError,IndexError,): # TypeError for int p, IndexError for numpy scalar (new scipy) self.parameters = [p] self.message = msg
[docs] def f_factory(self): """Stub for fit function factory, which returns the fit function. Override for derived classes. """ def fitfunc(p,x): # return f(p,x); should be a numpy ufunc raise NotImplementedError("base class must be extended for each fit function") return fitfunc
def _get_initial_values(self, parameters=None): p0 = numpy.asarray(self.initial_values()) if parameters is not None: try: p0[:] = parameters except ValueError: raise ValueError("Wrong number of custom initital values {0!r}, should be something like {1!r}".format(parameters, p0)) return p0
[docs] def initial_values(self): """List of initital guesses for all parameters p[]""" # return [1.0, 2.0, 0.5] raise NotImplementedError("base class must be extended for each fit function")
[docs] def fit(self,x): """Applies the fit to all *x* values""" fitfunc = self.f_factory() return fitfunc(self.parameters,numpy.asarray(x))
[docs]class FitExp(FitFunc): """y = f(x) = exp(-p[0]*x)""" def f_factory(self): def fitfunc(p,x): return numpy.exp(-p[0]*x) # exp(-B*x) return fitfunc def initial_values(self): return [1.0] def __repr__(self): return "<FitExp "+str(self.parameters)+">"
[docs]class FitExp2(FitFunc): """y = f(x) = p[0]*exp(-p[1]*x) + (1-p[0])*exp(-p[2]*x)""" def f_factory(self): def fitfunc(p,x): return p[0]*numpy.exp(-p[1]*x) + (1-p[0])*numpy.exp(-p[2]*x) return fitfunc def initial_values(self): return [0.5,0.1,1e-4] def __repr__(self): return "<FitExp2"+str(self.parameters)+">"
[docs]class FitLin(FitFunc): """y = f(x) = p[0]*x + p[1]""" def f_factory(self): def fitfunc(p,x): return p[0]*x + p[1] return fitfunc def initial_values(self): return [1.0,0.0] def __repr__(self): return "<FitLin"+str(self.parameters)+">"
[docs]class FitGauss(FitFunc): r"""y = f(x) = p[2] * 1/sqrt(2*pi*p[1]**2) * exp(-(x-p[0])**2/(2*p[1]**2)) Fits an un-normalized gaussian (height scaled with parameter p[2]). * p[0] == mean $\mu$ * p[1] == standard deviation $\sigma$ * p[2] == scale $a$ """ def f_factory(self): def fitfunc(p,x): return p[2] * 1.0/(p[1]*numpy.sqrt(2*numpy.pi)) * numpy.exp(-(x-p[0])**2/(2*p[1]**2)) return fitfunc def initial_values(self): return [0.0,1.0,0.0] def __repr__(self): return "<FitGauss"+str(self.parameters)+">"