Krioma.net

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:

#!/usr/bin/python #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
     fitfun = lambda p, x: p[0] + p[1] * x # define fitting function
     errfun = lambda p, x, y, yerr: (y - fitfun(p, x)) / yerr
     p0 = [0.0, 0.0]
     out = optimize.leastsq(errfun, p0, args=(logx, logy, logyerr), full_output=1)
     index = out[0][1]
     amp = 10.0**out[0][0]
     index_Err = sqrt( out[1][0][0] )
     amp_Err = sqrt( out[1][1][1] ) * amp
     r = spearcor(logx,logy) #Spearman rank correlation test
     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]