#!/usr/bin/env python3
#
'''
.. module:: expfit
:synopsis: Fit exponential function a*exp*b*x+c
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import numpy,scipy.optimize
def expfunc(p,x):
f=p[0]*numpy.exp(p[1]*x)+p[2]
return f
[docs]def expfit(x,y,yerr=None,debug=False):
'''
fit exponential function y=a*exp(b*x)+c
Args:
* `x` = x values
* `y` = y values to be fitted
Kwargs:
* `yerr` = sigma errors of the y-values
* `debug` = when True, print/plot debug information when fit does not converge
Returns:
* `p` = fitted parameter values for the function y=p[0]*exp(p[1]*x)+p[2]
* `yfit` = fitted function values
'''
hfunc=lambda x, p1, p2, p3: expfunc((p1,p2,p3),x)
# starting values for fit
ih=min((x.size-1),(x.size-x.size//4))
p=numpy.zeros(3,dtype=float)
rc=(y[-1]-y[0])/(x[-1]-x[0])
p[1]=1.0/y[ih]
p[0]=rc/2.7/p[1]
p[2]=-p[0]
if yerr is None:
yerr=numpy.ones(y.size,dtype=float)
try:
pf,covar = scipy.optimize.curve_fit(hfunc,x,y,p.copy(),yerr)
except RuntimeError:
print("scipy.optimize RuntimeError in function expfit")
print("x: ",x)
print("y: ",y)
print("initial values: ",p)
ih=max((0,(ih-1)))
p[1]=1.0/y[ih]
p[0]=1.0/3.0/p[1]
p[2]=-p[0]
try:
pf,covar = scipy.optimize.curve_fit(hfunc,x,y,p.copy(),yerr)
except RuntimeError:
print("scipy.optimize RuntimeError in function expfit")
print("initial values: ",p)
if debug:
import matplotlib.pyplot as plt
f,ax=plt.subplots(1)
ax.plot(x,y,'bo')
yfit=expfunc(p,x)
ss=x.argsort()
ax.plot(x[ss],yfit[ss],'r-')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Exponential fit attempt failed')
plt.show()
plt.close('all')
return (None,None)
yfit=expfunc(pf,x)
return (pf,yfit)
# ===============================================================
if __name__ == '__main__':
import matplotlib.pyplot as plt
y=numpy.array([1486.,1739.,5412.,5895.,5945.,6489.,8041.,8905.])
x=numpy.array([5097.,6109.,19507.,21085.,21253.,22967.,27615.,30020.])
p,yfit=expfit(x,y)
xx=numpy.arange(int(1.25*x.max()))
yy=expfunc(p,xx)
f,ax=plt.subplots(2)
ax[0].plot(x,y,'bo',xx,yy,'r-')
ax[0].set_title('Test of fitted expfunc')
ax[1].plot(x,y-yfit,'bo')
plt.show()
plt.close('all')