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

user_13516182 0
Viewed 781

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!

best-fit-curve lmfit model-fitting scipy-optimize curve-fitting
3 Answers

I've got two ideas, and I think the first one is probably your culprit.

  1. Using nan_policy = 'omit' seems to me to work in very specific cases. If you get the error message "nan_policy = 'omit' will probably not work" when trying to fit, well, it probably won't work. You could do a simple check for NaN values to confirm that the function outputs NaN values for your interval.
  2. The bounds for the variables are huge. Try raising the minimum for the intervals.
user_129965390

Without real data and a complete example, it is very hard to guess what might be going wrong. So, this will include some advice on how to approach the problem.

First: since you are doing curve-fitting, and have a model function, I recommend starting with your 2nd version, using lmfit.Model. But, I would suggest explicitly making a set of Parameters, as with:

from lmfit import Model
def A(t, gamma, Rp, Rs, a, orb_inclination, Rin, Rout):  
    Alpha = np.zeros(len(t))
    Xp, Zp = planetary_position(t, a, orb_inclination)
    values_rho = rho(Xp, Zp, Rs)
    v_W11 = W11(Xp,Zp, gamma, Rs)
    v_W11 = pd.Series(v_W11)
    v_zeta = zeta(Xp,Zp, gamma,Rs)
    v_zo = zo(Xp, Zp, gamma,Rs)
    v_W3 = W3(Xp,Zp, gamma,Rs)
    for i in range(len(values_rho)):
        Alpha[i]=np.where(values_rho[i]<1-gamma, np.pi*gamma**2*(1-u1-u2*(2-values_rho[i]**2-gamma**2/2)+(u1+2*u2)*v_W11[i] ) ,  np.where(((1-gamma<=values_rho[i]) and (values_rho[i]<=1+gamma)),  (1-u1-3*u2/2)*(gamma**2*np.arccos(v_zeta[i]/gamma)+np.arcsin(v_zo[i])-values_rho[i]*v_zo[i])+(u2/2)*values_rho[i]*((values_rho[i]+2*v_zeta[i])*gamma**2*np.arccos(v_zeta[i]/gamma) -v_zo[i]*(values_rho[i]*v_zeta[i]+2*gamma**2))+(u1+2*u2)*v_W3[i]    , 0)) 
    return Alpha

model = Model(model)
params = model.make_params(Rp=4.5*10**9, Rs=4.3*10**10, 
                           a=1.4*10**13, orb_inclination=1.2, 
                           Rin=4.5*10**9, Rout=13.5*10**9, tau=1)
result = gmodel.fit(y, params, t=t)
print(result.fit_report())

By itself, this won't solve the problem, but clarity counts. But, you can call that model function yourself or do 

  gmodel.eval(params, t=t)

and see what it actually calculates for any set of parameter values.

Second: you should be cautious about have variables in a fitting problem that span many orders of magnitude. Have the variables be more like order 1 (or, well between order 1.e-6 and 1.e6), and then multiply by factors of 1e9 or 1e12 as appropriate - or just work in units with values closer to 1. The numerics of fitting are all in double precision floating point, and the relative values of the parameters matters.

Third: your model function, yikes. Readability count. Writing an incomprehensible function does not help anyone. Including you. I guarantee that you do not know what this does. For example, you might be able to avoid the loop and just use the ufunc-ness of numpy, but it is impossible to tell. And to be clear, it is impossible to tell because you wrote it this way. Like what the heck are u1 and u2 supposed to be? Really, this function did not exist and you wrote a complete mess and then something went wrong.

So: write your model function as if you expect to read it next year, and then test what it calculates with reasonable input values. When that works, the fit should work too.

user_51797480

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.

user_135161820
Related