Gaussian Fitting in python

I spend a lot of my time working on noise statistics and of course and an important part of this is how to fit signals. I was asked earlier for an example code on how to fit a Gaussian, in particular fitting well defined signals. In python this is quite straight forward as you already have a bunch of nice libraries to count on.
So the first thing todo is to ensure you have these libraries. These are scipy, numpy, matplotlib and asciidata. Of course, you don’t need matplotlib if you don’t want todo the plotting but or asciddata if you don’t want to import an ascii file. (note: to students in 3rd year obs lab at Brum – you have this on hydra)
If you want to create some test data you can do something like:

from pylab import *
gaussian = lambda x: 1*exp(-(3-x)**2/10.) #change the parameters as you see fit
data = gaussian(arange(100))
x_pos = arange(100)

Now on to the fitting, this is simplistic and works well with “clean” data.

import warnings
warnings.simplefilter(“ignore”,DeprecationWarning)
import matplotlib
matplotlib.use(‘Agg’) #avoid X11 issues.
import scipy, numpy, pylab, asciidata
from numpy import *
from pylab import *
from scipy import *
from scipy.optimize import leastsq
gauss_fit = lambda p, x: p[0]*(1/sqrt(2*pi*(p[2]**2)))*exp(-(x-p[1])**2/(2*p[2]**2)) #1d Gaussian func
e_gauss_fit = lambda p, x, y: (gauss_fit(p,x) -y) #1d Gaussian fit
data = “test_gauss.dat” #your input data – should be x,y
demo = asciidata.open(data)
x_pos = double(demo[0]) #this is expected to be pixel number
y_power = double(demo[1])
v0= [1.0,mean(x_pos),1.0] #inital guesses for Gaussian Fit. $just do it around the peak
out = leastsq(e_gauss_fit, v0[:], args=(x_pos, y_power), maxfev=100000, full_output=1) #Gauss Fit
v = out[0] #fit parammeters out
covar = out[1] #covariance matrix output
xxx = arange(min(x_pos),max(x_pos),x_pos[1]-x_pos[0])
ccc = gauss_fit(v,xxx) # this will only work if the units are pixel and not wavelength
fig = figure(figsize=(9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(x_pos,y_power,’gs’) #spectrum
ax1.plot(xxx,ccc,’b–‘) #fitted spectrum
ax1.axvline(x=xxx[where(ccc == max(ccc))[0]][0],color=’r’) #max position in data
setp(gca(), ylabel=”power”, xlabel=”pixel position”)
pylab.savefig(“plotfitting.png”)
print “p[0], a: “, v[0]
print “p[1], mu: “, v[1]
print “p[2], sigma: “, v[2]

Here is the output plot (red line is the maximum, blue is the fit and green points are the data):
plotfitting
This code can all be found in a more useful format [simplegaussfit.py]
Download this, change the input data parameter to your asciifile and run as simplegaussfit.py – ensuring you have the correct file permissions to execute (you can change this with chmod).
Related to this page is another page I had on [fitting arbitrary functions to data].

Advertisements

2 responses to “Gaussian Fitting in python”

  1. Krioma.net Blog says :

    Multiple Gaussian Fitting in Python

    Yesterday I showed you [how to fit a single Gaussian in some data]. Today lets deal with the case of two Gaussians. This came about due to some students trying to fit two Gaussian’s to a shell star as the…

  2. Krioma.net Blog says :

    Spectral Extraction in Python

    I’ve wrote about how to read in FITS files in Python before, but I thought I’d readdress as I’ve been writing lots about fitting and wanted to build up to fitting properly calibrated data. So in this example I will…

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: