Archive | February 2011

Fitting data with python

As an astronomer I do lots of data fitting with python (utilising scipy) and thought it was about time I’d put up one of my basic procedures:

#by Samuel George; 03/02/2011
import matplotlib; matplotlib.use(‘Agg’) #avoid X11 issues.
import scipy, pylab, asciidata
from pylab import *
from scipy import optimize
from scipy.stats.stats import spearmanr as spearcor
def fitting_func(logx, logy, logyerr): # fitting
&nbsp&nbsp&nbsp&nbsp fitfun = lambda p, x: p[0] + p[1] * x # define fitting function
&nbsp&nbsp&nbsp&nbsp errfun = lambda p, x, y, yerr: (y – fitfun(p, x)) / yerr
&nbsp&nbsp&nbsp&nbsp p0 = [0.0, 0.0]
&nbsp&nbsp&nbsp&nbsp out = optimize.leastsq(errfun, p0, args=(logx, logy, logyerr), full_output=1)
&nbsp&nbsp&nbsp&nbsp index = out[0][1]
&nbsp&nbsp&nbsp&nbsp amp = 10.0**out[0][0]
&nbsp&nbsp&nbsp&nbsp index_Err = sqrt( out[1][0][0] )
&nbsp&nbsp&nbsp&nbsp amp_Err = sqrt( out[1][1][1] ) * amp
&nbsp&nbsp&nbsp&nbsp r = spearcor(logx,logy) #Spearman rank correlation test
&nbsp&nbsp&nbsp&nbsp return out[0],out[1],index,amp,index_Err,amp_Err,r #p[0], p[1], covariance matrix, index of powerlaw, amp powerlaw, errors+spearman rank
pout, covar,index,amp,indexErr,ampErr,r = fitting_func(log10(xdata), log10(ydata), (errory / log10(ydata)))

This code can also be found [here]