[SOLVED] lmfit and scipy curve_fit return initial guesses as best-fitted parameters

lmfit and scipy curve_fit return initial guesses as best-fitted parameters

I want to fit a function to some data and I’ m facing a problem. I’ ve tried to use lmfit or curve_fit from scipy. Below I describe the problem.

Here is my data:

dataOT = pd.read_csv("KIC3239945e.csv", sep=';') 
t=dataOT['time']
y=dataOT['flux']

Also, here is the model-function to be fitted to the data:

def model(t, Rp, Rs, a, orb_inclination, Rin, Rout, tau):  
    gps=Rp/Rs
    gis=Rin/Rs
    gos=Rout/Rs
    Agps=A(t, gps, Rp, Rs, a, orb_inclination, Rin, Rout)
    Agos=A(t, gos, Rp, Rs, a, orb_inclination, Rin, Rout)
    Agis=A(t, gis, Rp, Rs, a, orb_inclination, Rin, Rout)
    return (np.pi*(1-u1/3-u2/6)-Agps-(1-np.exp(-tau))*(Agos-Agis))/(np.pi*(1-u1/3-u2/6))

where u1, u2 are known numbers and the parameters to be fitted are: Rp, Rs, a, orb_inclination, Rin, Rout, tau and they are contained in the quantities Agps, Agos, Agis. Here is the definition of function A:

def A(t, gamma, Rp, Rs, a, orb_inclination, Rin, Rout):  
Xp,Zp= planetary_position(t, a, orb_inclination)
return np.where(rho(Xp,Zp,Rs)<1-gamma,
                np.pi*gamma**2*(1-u1-u2*(2-rho(Xp,Zp,Rs)**2-gamma**2/2)+(u1+2*u2)*W11(Xp,Zp,gamma,Rs) ) , 
                np.where(np.logical_and(  (1-gamma<=rho(Xp,Zp,Rs)),  (rho(Xp,Zp,Rs)<=1+gamma)  ), 
                (1-u1-3*u2/2)*(gamma**2*np.arccos(zeta(Xp,Zp,gamma,Rs)/gamma)+np.arcsin(zo(Xp,Zp,gamma,Rs))-rho(Xp,Zp,Rs)*zo(Xp,Zp,gamma,Rs))+(u2/2)*rho(Xp,Zp,Rs)*((rho(Xp,Zp,Rs)+2*zeta(Xp,Zp,gamma,Rs))*gamma**2*np.arccos(zeta(Xp,Zp,gamma,Rs)/gamma)-zo(Xp,Zp,gamma,Rs)*(rho(Xp,Zp,Rs)*zeta(Xp,Zp,gamma,Rs)+2*gamma**2))  +(u1+2*u2)*W3(Xp,Zp,gamma,Rs)    , 0))

1st attempt: curve_fit

from scipy.optimize import curve_fit
p0=[4.5*10**9, 4.3*10**10, 1.4*10**13, 1.2, 4.5*10**9, 13.5*10**9, 1]
popt, pcov = curve_fit(model, t, y, p0, bounds=((0, 0, 0, 0, 0, 0 ,0 ),(np.inf, np.inf, np.inf,np.inf, np.inf, np.inf ,np.inf )), maxfev=6000)
print(popt)

2nd attempt: lmfit

   from lmfit import Parameters, Minimizer, report_fit, Model
gmodel=Model(model)

def residual(p,t, y):
    Rp=p['Rp']
    Rs=p['Rs']
    a=p['a']
    orb_inclination=p['orb_inclination']
    Rin=p['Rin']
    Rout=p['Rout']
    tau=p['tau']
    tmp = model(t, Rp, Rs, a, orb_inclination, Rin, Rout, tau)-y
    return tmp

p = Parameters()

p.add('Rp' ,  value=0.000394786,     min= 0,max=1)
p.add('Rs' ,  value=0.003221125,    min= 0,max=1)
p.add('a',   value=1.86,            min= 0,max= 1)
p.add('orb_inclination',  value=1,   min= 0,max=4)
p.add('Rin',  value=0.000394786,    min= 0,max=1)
p.add('Rout',  value=0.000394786,    min= 0,max=1)
p.add('tau',  value=0,                 min= 0,max=2)

mini = Minimizer(residual,params=p,fcn_args=(t,y))

out = mini.minimize(method='leastsq')

print(report_fit(out))

All cases return as best-fitted parameters the initial guesses. What should I do in order to make it work properly?

Note:Assuming that the parameters are known the model has the expected behavior (Figure 1), so I suppose that the model is well-defined and the problem is not related with this.

Any help would be appreciated. Thank you in advance!

Solution

I solved the problem by reducing the number of parameters. Also, another problem was that one of the parameters was not affecting the fitting at all.