Krioma.net

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