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 spectral line was altered from a simple Gaussian, actually there is a nice P-Cygni dip in there data so you should be able to recover the absorption line by this kind of fitting. Anyway, fitting 2 Gaussian’s is basically the same thing as fitting one in python but with the added function. I’m thinking I’ll work up how to deal with a whole spectrum including source finding. In this case what you have to deal with is that there are two sources and so a rough estimation of the peak position of both is crucial to the fit (well in the way it is implemented).
You basically want to end up with something like this:
Fit 2 Gaussians in Python
Where the green points are the data, the blue dashed line the fit and the red line where the maximum is in the array. The fitting is done via a least-squares fitting routine.
In this example we will first start by generating some data, skipping the input from a file (but you can use the code from the 1 Gaussian example).

import scipy, numpy, pylab, asciidata
from numpy import *
from pylab import *
from scipy import *
from scipy.optimize import leastsq
#generate some data
gaussian = lambda x: 3*exp(-(10-x)**2/10.) + 1*exp(-(30-x)**2/10.)#change the parameters as you see fit
y_power = gaussian(arange(100))
x_pos = arange(100)
gauss_fit = lambda p, x: p[0]*(1/sqrt(2*pi*(p[2]**2)))*exp(-(x-p[1])**2/(2*p[2]**2))+p[3]*(1/sqrt(2*pi*(p[5]**2)))*exp(-(x-p[4])**2/(2*p[5]**2)) #1d Gaussian func
e_gauss_fit = lambda p, x, y: (gauss_fit(p,x) -y) #1d Gaussian fit
v0= [1,10,1,1,30,1] #inital guesses for Gaussian Fit. – just do it around the peaks
out = leastsq(e_gauss_fit, v0[:], args=(x_pos, y_power), maxfev=100000, full_output=1) #Gauss Fit
v = out[0] #fit parameters 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”)
print “p[0], a1: “, v[0]
print “p[1], mu1: “, v[1]
print “p[2], sigma1: “, v[2]
print “p[3], a2: “, v[0]
print “p[4], mu2: “, v[1]
print “p[5], sigma2: “, v[2]

The full detailed listing is given in []

One response to “Multiple Gaussian Fitting in Python”

  1. 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: Logo

You are commenting using your 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: