curve fitting - Not able to fit a function with scipy.optimize.curve_fit() -
i fit following function:
def invlaplace_stehfest2(time,el,tau): upsilon=0.25 pmax=6.9 e0=0.0154 m=8 results=[] t in time: func=0 k in range(1,2*m+1): sum=0 j in range(int(math.floor((k+1)/2)),min(k,m)+1): dummy=j**(m+1)*scipy.special.binom(m,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(m) sum+=dummy s=k*math.log(2)/t[enter image description here][1] func+=(-1)**(m+k)*sum*pmax*el/(mp.exp(tau*s)*mp.expint(1,tau*s)*e0+el)/s func=func*math.log(2)/t results.append(func) return [float(i) in results] to use following data:
data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) with folowing guess:
guess=np.array([0.1,0.05]) and scipy.optimize.curve_fit() folow:
parameter,covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) for reason don't understand not able fit data correctly. following graph.
my function undoubtedly correct since when use correct guess:
guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) i able fit correctly data.
any hint on why not able fit correctly? should use method?
here hole program:
############################################## import matplotlib import matplotlib.pyplot plt import matplotlib.mlab mlab import matplotlib.pylab mplab import math math import * import numpy np import scipy scipy.optimize import curve_fit import mpmath mp ############################################################################################ def invlaplace_stehfest2(time,el,tau): upsilon=0.25 pmax=6.9 e0=0.0154 m=8 results=[] t in time: func=0 k in range(1,2*m+1): sum=0 j in range(int(math.floor((k+1)/2)),min(k,m)+1): dummy=j**(m+1)*scipy.special.binom(m,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(m) sum+=dummy s=k*math.log(2)/t func+=(-1)**(m+k)*sum*pmax*el/(mp.exp(tau*s)*mp.expint(1,tau*s)*e0+el)/s func=func*math.log(2)/t results.append(func) return [float(i) in results] ############################################################################################ ###constant### ####value#### data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) ###fitting### guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) #guess=np.array([0.1,0.05]) parameter,covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) print parameter residu=sum(data_relax-invlaplace_stehfest2(data_time,parameter[0],parameter[1])) graph_curves=plt.figure() ax = graph_curves.add_subplot(111) ax.plot(data_time,invlaplace_stehfest2(data_time,parameter[0],parameter[1]),"-") ax.plot(data_time,data_relax,"o") plt.show()
non-linear fitters such default levenberg-marquardt solver used in scipy.optimize.curve_fit(), iterative solvers, can stop in local minimum in error space. if error space "bumpy" initial parameter estimates become important, in case.
scipy has added differential evolution genetic algorithm optimize module, can used determine initial parameter estimates curve_fit(). scipy's implementation uses latin hypercube algorithm ensure thorough search of parameter space, requiring parameter upper , lower bounds within search. can see below, have used scipy module replace hard-coded values value named "guess" in code. not run quickly, slower correct result better fast wrong result. try code:
import matplotlib import matplotlib.pyplot plt import matplotlib.mlab mlab import matplotlib.pylab mplab import math math import * import numpy np import scipy scipy.optimize import curve_fit scipy.optimize import differential_evolution import mpmath mp ############################################################################################ def invlaplace_stehfest2(time,el,tau): upsilon=0.25 pmax=6.9 e0=0.0154 m=8 results=[] t in time: func=0 k in range(1,2*m+1): sum=0 j in range(int(math.floor((k+1)/2)),min(k,m)+1): dummy=j**(m+1)*scipy.special.binom(m,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(m) sum+=dummy s=k*math.log(2)/t func+=(-1)**(m+k)*sum*pmax*el/(mp.exp(tau*s)*mp.expint(1,tau*s)*e0+el)/s func=func*math.log(2)/t results.append(func) return [float(i) in results] ############################################################################################ ###constant### ####value#### data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) # function genetic algorithm minimize (sum of squared error) def sumofsquarederror(parametertuple): return np.sum((data_relax - invlaplace_stehfest2(data_time, *parametertuple)) ** 2) ###fitting### #guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) #guess=np.array([0.1,0.05]) parameterbounds = [[0.0, 1.0], [0.0, 10.0]] # "seed" numpy random number generator repeatable results # note ".x" here return parameter estimates in example guess = differential_evolution(sumofsquarederror, parameterbounds, seed=3).x parameter,covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) print parameter residu=sum(data_relax-invlaplace_stehfest2(data_time,parameter[0],parameter[1])) graph_curves=plt.figure() ax = graph_curves.add_subplot(111) ax.plot(data_time,invlaplace_stehfest2(data_time,parameter[0],parameter[1]),"-") ax.plot(data_time,data_relax,"o") plt.show()
Comments
Post a Comment