# 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):

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].

Multiple Gaussian Fitting in PythonYesterday 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…

Spectral Extraction in PythonI’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…