Search code examples
pythongaussiandata-fitting

Fitting data with multiple Gaussian profiles in Python

I have some data (data.txt) and am trying to write a code in Python to fit them with Gaussian profiles in different ways to obtain and compare the peak separation and the under curve area in each case:

  1. with two Gaussian profiles (considering the little peaks on top and ignoring the shoulders; the red profiles)
  2. with two Gaussian profiles (ignoring the little peaks on top and considering the whole single peak at top and the shoulders; the black profiles)
  3. with three Gaussian profiles (considering a tall peak on the two shorter ones in the shoulders; the green profiles)

I tried several scripts, but I failed in all.

The profiles in these plots are fake, and I just added them to illustrate better what I mean.

plots

Solution

  • One approach to this is as follows:

    1. Define the function you want to fit to the data, i.e. the sum of all components that should be in there. In your case this is multiple Gaussians.
    2. Find initial guesses for your parameters.
    3. Fit your fitting function to the data, using a strategy to your liking.

    I took a go at your data, and below is a very simple example of fitting for three Gaussian components and a continuum offset, using SciPy's curve_fit method. I'll leave the rest to you. This should allow you to figure out the other cases as well. Note that initial guesses generally are important, so it's best to make an educated guess somehow, to get as close to the optimal value as possible.

    Code

    from scipy.optimize import curve_fit

import matplotlib.pyplot as plt
import numpy as np

def gaussian(x, A, x0, sig):
    return A*np.exp(-(x-x0)**2/(2*sig**2))

def multi_gaussian(x, *pars):
    offset = pars[-1]
    g1 = gaussian(x, pars[0], pars[1], pars[2])
    g2 = gaussian(x, pars[3], pars[4], pars[5])
    g3 = gaussian(x, pars[6], pars[7], pars[8])
    return g1 + g2 + g3 + offset

vel, flux = np.loadtxt('data.txt', unpack=True)
# Initial guesses for the parameters to fit:
# 3 amplitudes, means and standard deviations plus a continuum offset.
guess = [4, -50, 10, 4, 50, 10, 7, 0, 50, 1]
popt, pcov = curve_fit(multi_gaussian, vel, flux, guess)

plt.figure()
plt.plot(vel, flux, '-', linewidth=4, label='Data')
plt.plot(vel, multi_gaussian(vel, *popt), 'r--', linewidth=2, label='Fit')
plt.legend()
plt.show()

    Result

    Three-Gaussian fit