#!/usr/bin/env python3
#
'''
.. module:: gaussfit
:synopsis: Fit guassian function to data
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import sys
import numpy
import scipy.optimize
[docs]def gauss(p,x):
'''
Compute Gauss function
Args:
* `p` = Gauss parameters a*exp(-(x-b)**2/(2*c**2))
* `p[0]` = a - scaling
* `p[1]` = b - position
* `p[2]` = c - sigma width
* `x` = x-axis values
Returns:
* `y` = function values
'''
y=p[0]*numpy.exp(-(x-p[1])**2/(2.0*p[2]**2))
return y
[docs]def gaussfit(x,y):
'''
Fit Gauss function to data.
use Chi-square statistics fit
Args:
* `x` = x-axis values
* `y` = data values to be fitted
Returns:
* `p` = list of fitted gauss parameters
* `p[0]` = a - scaling
* `p[1]` = b - position
* `p[2]` = c - sigma width
* `perr` = list of error values for the fitted parameters
* `perr[0]` = a - scaling error
* `perr[1]` = b - position error
* `perr[2]` = c - sigma width error
* `fit` = fitted function values
'''
gaussfunc=lambda x, p1, p2, p3: gauss((p1,p2,p3),x)
p=numpy.zeros(3)
# initial fit values
box=numpy.array([1,2,3,4,3,2,1],dtype=float)
box=box/numpy.sum(box)
ccf=numpy.convolve(y,box,mode='same')
mm=numpy.argmax(y)
p[0]=ccf[mm]
p[1]=x[mm]
indx,=numpy.where(ccf > (0.5*p[0]))
if indx.size > 1:
p[2]=x[indx[-1]]-x[indx[0]]
else:
p[2]=x[1]-x[0]
yy=numpy.where(y <= 0.0,1.0,y)
ye=numpy.sqrt(yy)
try:
p1, covar =scipy.optimize.curve_fit(gaussfunc, x, y, p.copy(), ye) # first fit
except:
print("Fit error:", sys.exc_info()[0])
print("Error, no proper fit for gaussfit")
return (p,numpy.zeros(len(p)),gauss(p,x))
ye=numpy.sqrt(gauss(p1,x))
# reiterate with errors based on model (Chi-square statistics)
try:
p2, covar =scipy.optimize.curve_fit(gaussfunc, x, y, p1.copy(), ye)
ye=numpy.sqrt(gauss(p1,x))
p2, covar =scipy.optimize.curve_fit(gaussfunc, x, y, p2.copy(), ye)
except:
print("C-stat fit error:", sys.exc_info()[0])
print("Cannot fit C-stat for gaussfit")
return (p1,numpy.zeros(len(p1)),gauss(p1,x))
perr = numpy.sqrt(numpy.diag(covar)) # parameter error estimates
return (p2,perr,gauss(p2,x))
#================================================================================
if __name__ == "__main__": # test program
import matplotlib
import matplotlib.pyplot as plt
x=numpy.arange(500)/100.0-2.5 # x-range -2.5 2.5
p=[10.0,0.0,0.1]
p2=[3.0,2.0,0.5]
y=gauss(p,x)+numpy.random.random(x.size)*2.0+gauss(p2,x)
pp,pe,fit=gaussfit(x.astype(float),y.astype(float))
print("Fitpars: ",pp)
print("Errors: ",pe)
f,ax=plt.subplots(1)
ax.plot(x,y,'b-')
ax.plot(x,fit,'r-')
ax.set_title('Gauss fit test')
plt.show()
plt.close('all')