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’re 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:
#defines model function for the data to be fitted to def gauss(x, *p): 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)