Gaussian Fit in Python

Created on 2014-11-30 (Sun)

Gaussian fitting in Python... what a pain. Have you ever wanted to fit some data to a Gaussian in Python in 2 minutes? I did. It's easy right... there must be a function that does least square fits... hopefully there's some instructions on how to use it. Call it with a Gaussian as the work function and you are done. Or so I thought. I also thought it would be smart to get on stackoverflow.com and see what code lurked around there. So this is the story of what I found:

I started with this incorrect - just plain wrong version (note to self 1: post a message correcting it). Then I got to this version which is confusing if you want to get a Gaussian in 2 minutes. Though this last version had an interesting side link - which even though not about Gaussians is somewhat relevant to fitting functions. Finally I get to this correct starting point - which was a non-fitted version of what I wanted (note to self 2: in the future use the Cookbook more than stackoverflow). Then once I got confident I was on the right path there's this bug related to using floats - a strange apparition waiting to ruin my day.

So without much further ado here is - finally - my correct and functional version:

def gauss(x, *p):
   '''defines model function for the data to be fitted to'''
   A, mu, sigma = p
   # I need to cast denominator and numerator as float64 because of a 
   # bug in numpy
   numer = np.float64(-(x-mu)**2.)
   denom = np.float64((2.*sigma**2.))
   retVal = A*np.exp(numer/denom)
   return retVal

def fitToGaussian(xData, yData):
   A = yData.max()
   mean = (sum(xData*yData)/sum(yData))
   sigma =((sum(yData*(xData-mean)**2)/sum(yData)))
   print ('A:' + str(A) + '; mean: ' + str(mean) + 
          '; sigma: ' + str(sigma))
   # p0 is the initial guess for the fitting coefficients (A, mu and 
   # sigma)
   p0 = [A, mean, sigma]
   coeff, var_matrix = curve_fit(gauss, xData, yData, p0=p0)
   return coeff

dataX = loadDataX()
dataY = loadDataY()

coeff = fitToGaussian(dataX, dataY)