0

I have a set of points in the first quadrant that look like a gaussian, and I am trying to fit it using a gaussian in python and my code is as follows:

import pylab as plb
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import asarray as ar,exp
import math


x=ar([37,69,157,238,274,319,391,495,533,626,1366,1855,2821,3615,4130,4374,6453,6863,7021,
    7951,8646,9656,10464,11400])
y=ar([1.77,1.67,1.65,1.17,1.34,1.46,0.75,1,0.8,1.02,0.65,0.69,0.44,0.44,0.55,0.43,0.75,0.27,0.26,
    0.44,0.04,0.44,0.26,0.04])



n = 24                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = math.sqrt(sum(y*(x-mean)**2)/n)        #note this correction

def gaus(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))

popt,pcov = curve_fit(gaus,x,y,p0=None, sigma=None)  #'''p0=[1,mean,sigma]'''

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gaus(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

And the output is: this figure: http://s2.postimg.org/wevggkc95/Workspace_1_022.png

Why are all the red points coming below, Also note that I am interested in a half gaussian as my data is like that, so my y values are big at first and then decreasing like one side of the gaussian bell. Can anyone tell me how to fit this curve in python, (in case it cannot be fit to gaussian). Or in other words, I want code to fit the half(left side) gaussian of my points (in the first quadrant only). Note that my points cannot be fit as an exponentially decreasing curve as I tried that earlier, and it is not fitting well at lower 'x' values.

PsJain
  • 3
  • 4

1 Answers1

0

Apparently your data do not fit well or easily to a Gaussian function. You use the default initial guesses for p0 = [1,1,1] which is so far away from any kind of optimal choice that curve_fit gives up before it gets started (check the values of popt=[1,1,1] and pcov=[inf, inf, inf]). You could try with better guesses (e.g. p0 = [2,0, 2000]), but on my system it won't converge: Optimal parameters not found: Number of calls to function has reached maxfev = 800.

To fit a "half-Gaussian", don't float the centre position x0 (just leave it equal to 0):

def gaus(x,a,sigma):
    return a*exp(-(x)**2/(2*sigma**2))

p0 = [1.2, 4000]
popt,pcov = curve_fit(gaus,x,y,p0=p0)

enter image description here

Unless you have a particular reason for wanting to fit a Gaussian, why not do a more robust linear least squares fit to a polynomial, e.g.:

pfit = np.polyfit(x, y, 3)
poly = np.poly1d(pfit)
xnx
  • 24,509
  • 11
  • 70
  • 109
  • Thanks, I know it does not look like a gaussian, but it certainly is like the right part of the Gaussian. I wanted to ask if there could be a way to fit only such data sets to half gaussians if possible. – PsJain Nov 26 '15 at 13:25
  • @PsJain Ah, OK - I've edited to give you the way to do this. – xnx Nov 26 '15 at 13:30
  • Now that worked very well! Thank you very very much for your quick and precise response! – PsJain Nov 26 '15 at 13:40