1

Firstly this is an assignment I've been set so I'm only after pointers, and I am restricted to using the following libraries, NumPy, SciPy and MatPlotLib.

We have been given a txt file which includes x and y data for a resonance experiment and have to fit both a gaussian and lorentzian fit. I'm working on the gaussian fit at the minute and have tried following the code laid out in a previous question as a basis for my own code. (Gaussian fit for Python)

from numpy import *
from matplotlib import *
import matplotlib.pyplot as plt  ##Import Libraries
import pylab
from scipy.optimize import curve_fit

####################################

energy,intensity=numpy.loadtxt('resonance_data.txt',unpack=True)
print energy
print intensity  ##Load in text file and print the arrays




#####################################

n = size(energy)
mean = 30.7
sigma = 10
intensity0 = 45

def gaus(energy,intensity0,energy0,sigma):
       return intensity0*exp(-(energy-energy0)**2/(sigma**2))

popt,pcov = curve_fit(gaus,energy,intensity,p0=[45,mean,sigma])



plt.plot(energy,intensity,'o')
plt.xlabel('Energy/eV')
plt.ylabel('Intensity')
plt.title('Plot of Intensity against Energy')  ##Plot raw data along with axis labels         and title
plt.plot(energy,gaus(energy,*popt))
plt.show()

Which returns the following graph

enter image description here

If I keep the expressions for mean and sigma, as in the url posted the curve fit is a horizontal line, so I'm guessing the problem lies in the curve fit not converging or something.

I am somewhat of a beginner so any pointers are helpful :)

Thanks

UPDATE: I've managed to improve the Gaussian fit and get the Lorentzian working, and even (I think) managed to calculate the sum of the residuals in each case.

Cheers guys!

Thanks again guys

Community
  • 1
  • 1
user3507401
  • 11
  • 1
  • 3
  • 1
    Somehow the pic you links presents itself as a decent fit.... – CT Zhu Apr 07 '14 at 17:52
  • The general shape is close to what I would expect but the height of the peak is quite far off. I'm aware its' pretty skewed so its not going to fit perfectly, and that's sort of the point of the assignment, but I would still expect better? – user3507401 Apr 07 '14 at 18:18
  • Are you supposed to be fitting a one or a two dimensional gaussian distribution? – Paul Apr 07 '14 at 18:57
  • I would use a lorentzian for this data. Think your fit is converged. Increasing the amplitude would probably decrease the fit of the wings. – tillsten Apr 07 '14 at 20:29

1 Answers1

0

Looks like your data skews heavily to the left, why Gaussian? Not Boltzmann, Log-Normal, or anything else?

Much of these are already implemented in scipy.stats. See scipy.stats.cauchy for lorentzian and scipy.stats.normal gaussian. An example:

import scipy.stats as ss
A=ss.norm.rvs(0, 5, size=(100)) #Generate a random variable of 100 elements, with expected mean=0, std=5
ss.norm.fit_loc_scale(A) #fit both the mean and std
(-0.13053732553697531, 5.163322485150271) #your number will vary.

And I think you don't need the intensity0 parameter, it is just going to be 1/sigma/srqt(2*pi), because the density function has to sum up to 1.

Gabriel
  • 37,092
  • 64
  • 207
  • 369
CT Zhu
  • 49,083
  • 16
  • 110
  • 125
  • Thanks for the reply, we've been told to plot both the Gaussian and Lorentzian and then perform statistical tests to show which is best, so that is the reason for the choice. I'll take a look at the scipy functions Cheers. – user3507401 Apr 07 '14 at 18:16
  • @user3507401 you can [have a look to this example here](http://stackoverflow.com/a/16651955/832621) to see how to test many different distributions using SciPy. – Saullo G. P. Castro Apr 08 '14 at 13:07