From b8af43314e05a2345d84cc248782e57f11674323 Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Sun, 2 Mar 2008 21:33:33 -0500 Subject: tweaks --- stats.py | 90 +++++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 58 insertions(+), 32 deletions(-) (limited to 'stats.py') diff --git a/stats.py b/stats.py index 9acf3f4..eea94a2 100644 --- a/stats.py +++ b/stats.py @@ -50,7 +50,7 @@ def rchisqr_poisson(tp,ffunc,tx,ty,floor=1.): terr = floored_root(tx,min=floor); return rchisqr(tp,ffunc,tx,ty,terr); -def fittodata(tx,ty,p0,ffunc,err=None,floor=1.,ret_all=False): +def fittodata(tx,ty,p0,ffunc,err=None,floor=1.,ret_all=False,param_names=None): """does a fit and returns a dict of shit. ret_all=True gives more shit.""" if err == None: @@ -101,50 +101,70 @@ def fittodata(tx,ty,p0,ffunc,err=None,floor=1.,ret_all=False): "err":err, "x":tx, "y":ty, - "success":success}; + "success":success, + "param_names":param_names}; else: return {"params":p1, "param_error":p1_err, - "rchisqr":rcs} + "rchisqr":rcs, + "param_names":param_names}; def fittodata_gaussian(tx,ty,p0,err=None,floor=1.,ret_all=False): """runs fittodata for a gaussian (first parameter is mean, - second is stddev, third is amplitude)""" + second is stddev)""" #gauss_func = lambda a,b: a[2]/(a[1] * sqrt(2*pi)) * exp(-.5 * ((b - a[0])/a[1])**2) #thing = std(ty); bigsum = scipy.integrate.simps(ty,x=tx) gauss_func = lambda a,b: bigsum/(a[1] * scipy.sqrt(2*pi)) * scipy.exp(-.5 * ((b - a[0])/a[1])**2) - return fittodata(tx,ty,p0,gauss_func,err=err,floor=floor,ret_all=ret_all); + return fittodata(tx,ty,p0,gauss_func,err=err,floor=floor,ret_all=ret_all, + param_names=["Mean","StdDev"]); + +def fittodata_gaussian_amp(tx,ty,p0,err=None,floor=1.,ret_all=False): + """runs fittodata for a gaussian with amplitude parameter (first parameter + is mean, second is stddev, third is amplitude)""" + bigsum = scipy.integrate.simps(ty,x=tx) + gauss_func = lambda a,b: a[2] * bigsum/(a[1] * scipy.sqrt(2*pi)) * scipy.exp(-.5 * ((b - a[0])/a[1])**2) + return fittodata(tx,ty,p0,gauss_func,err=err,floor=floor,ret_all=ret_all, + param_names=["Mean","StdDev","Amplitude"]); + +def fittodata_gaussian_plus(tx,ty,p0,err=None,floor=1.,ret_all=False): + """runs fittodata for a gaussian with extra parameters (first parameter + is mean, second is stddev, third is amplitude, fourth is constant offset, + fifth is linear trend)""" + bigsum = scipy.integrate.simps(ty,x=tx) + gauss_func = lambda a,b: a[2] * bigsum/(a[1] * scipy.sqrt(2*pi)) * scipy.exp(-.5 * ((b - a[0])/a[1])**2) + a[3] + a[4]*b + return fittodata(tx,ty,p0,gauss_func,err=err,floor=floor,ret_all=ret_all, + param_names=["Mean","StdDev","Amplitude","VertOffset","Slope"]); def fittodata_linear(tx,ty,p0,err=None,floor=1.,ret_all=False): """runs fittodata for a linear curve (first param is slope, second is vertical offset)""" lin_func = lambda a,b: a[0] * b + a[1] - return fittodata(tx,ty,p0,lin_func,err=err,floor=floor,ret_all=ret_all); + return fittodata(tx,ty,p0,lin_func,err=err,floor=floor,ret_all=ret_all, + param_names=["VertOffset (m)","Slope (b)"]); + +def fittodata_quadratic(tx,ty,p0,err=None,floor=1.,ret_all=False): + """runs fittodata for a quadratic curve (first param is a, + second is b, third is c: ax^2+bx+c)""" + quad_func= lambda a,b: a[0] * b*b + a[1]*b + a[2] + return fittodata(tx,ty,p0,quad_func,err=err,floor=floor,ret_all=ret_all, + param_names=["QuadaticTerm (a)","VertOffset (b)","Slope (c)"]); def fittodata_poisson(tx,ty,p0,err=None,floor=.5,ret_all=False): - """runs fittodata for a poisson (first param is the mean""" + """runs fittodata for a poisson (first param is the mean)""" bigsum = scipy.integrate.simps(ty,x=tx) if type(p0) == type(0.1): p0 = [p0,] poissfunc = lambda a,b: bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0]) - #def poissfunc(a,b): - # print "--------" - # print a - # print b - # print bigsum - # print scipy.power(a[0],b) - # print scipy.exp(a[0]) - # print bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0]) - # return bigsum * scipy.power(a[0],b) / scipy.factorial(b) * scipy.exp(-1 * a[0]) - return fittodata(tx,ty,p0,poissfunc,err=err,floor=floor,ret_all=ret_all); + return fittodata(tx,ty,p0,poissfunc,err=err,floor=floor,ret_all=ret_all, + param_names=["Mean"]); def fittodata_exp(tx,ty,p0,err=None,floor=1.,ret_all=False): """runs fittodata for an exponential (first param is x0, second param is tau, third param is vertical offset)""" - #exp_func = lambda a,b: a[0] * scipy.exp(- b/a[1]) + a[2] exp_func = lambda a,b: a[0] * e ** (- b/a[1]) + a[2] - return fittodata(tx,ty,p0,exp_func,err=err,floor=floor,ret_all=ret_all); + return fittodata(tx,ty,p0,exp_func,err=err,floor=floor,ret_all=ret_all, + param_names=["x0","tau","VertOffset"]); def fittodata_lorentz(tx,ty,p0,err=None,floor=1.,ret_all=False): """runs fittodata for a lorentzian (first param is gamma, @@ -153,7 +173,19 @@ def fittodata_lorentz(tx,ty,p0,err=None,floor=1.,ret_all=False): if type(p0) == type(0.1): p0 = [p0,] lorentz_func = lambda a,b: bigsum * a[0]/(2.*pi*((a[1]-b)**2+(a[0]/2.)**2. )) - return fittodata(tx,ty,p0,lorentz_func,err=err,floor=floor,ret_all=ret_all); + return fittodata(tx,ty,p0,lorentz_func,err=err,floor=floor,ret_all=ret_all, + param_names=["Gamma","Mean"]); + +def fittodata_lorentz_plus(tx,ty,p0,err=None,floor=1.,ret_all=False): + """runs fittodata for a lorentzian with extra parameters (first param + is gamma, second param is mean, third is constant offset, fourth is + linear trend, fifth is amplitute)""" + bigsum = scipy.integrate.simps(ty,x=tx) + if type(p0) == type(0.1): + p0 = [p0,] + lorentz_func = lambda a,b,: a[4]*bigsum * a[0]/(2.*pi*((a[1]-b)**2+(a[0]/2.)**2. )) + a[2] + a[3]*b + return fittodata(tx,ty,p0,lorentz_func,err=err,floor=floor,ret_all=ret_all, + param_names=["Gamma","Mean","VertOffset","Slope","Amplitude"]); def plotfit(indict,newfig=True): if not indict.has_key('x'): @@ -164,6 +196,7 @@ def plotfit(indict,newfig=True): y = indict['y'] err = indict['err'] ffunc = indict['ffunc'] + param_names = indict['param_names'] if newfig: pylab.figure(); pylab.errorbar(x,y,fmt='b.',yerr=err) d=pylab.plot(x,y,'k.') @@ -179,19 +212,12 @@ def plotfit(indict,newfig=True): p1_err = indict['param_error'] ax = pylab.gca() - extra_txt = 'reduced chi squared: %g'\ - % (indict['rchisqr']) + extra_txt = 'reduced chi^2: %g'% (indict['rchisqr']) for i in range(scipy.size(p1)): - extra_txt += '\np[%i]: %g +/- %g' \ - % (i,p1[i],p1_err[i]); - -# if(pylab.size(p1) == 2): -# extra_txt = 'reduced chi squared: %g\np[0]: %g\np[1]: %g' \ -# % (indict['rchisqr'],p1[0],p1[1]) -# else: -# extra_txt= 'reduced chi squared:%g\np[0]:%g, +/- %g\n p[1]:%g\np[2]:%g' \ -# % (indict['rchisqr'],p1[0],p1[1],p1[2]) - + if len(param_names)>i: + extra_txt += '\n%s: %.3g +/- %.2g' % (param_names[i],p1[i],p1_err[i]); + else: + extra_txt += '\np[%i]: %.3g +/- %.2g' % (i,p1[i],p1_err[i]); t = pylab.text(.77, .75, extra_txt, horizontalalignment='center', verticalalignment='center', -- cgit v1.2.3