2

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!

kop
  • 21
  • 5

3 Answers3

0

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.
sbjartmar
  • 23
  • 6
  • I changed the bounds of the variables. For example I wrote: `p.add('Rp' , value=4.5*10**9, min= 3*10**9,max=8*10**9)` but the same problem still exists. Also, neither residual function nor model function contain nan values. What else could be responsible for the wrong fitting-results? – kop May 11 '20 at 20:25
0

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.

M Newville
  • 7,486
  • 2
  • 16
  • 29
  • I changed the values of the parameters and so the bounds. Now, the 1st case I wrote before includes also: `params=gmodel.make_params(Rp=0.000394786, Rs=0.003221125, a=1.86, orb_inclination=1.2, Rin=0.000394786, Rout=0.000394786, tau=1)` `model_values=gmodel.eval(params, t=t)` When I run this, the following error appears: `ValueError: The model function generated NaN values and the fit aborted! Please check your model function and/or set boundaries on parameters where applicable. In cases like this, using "nan_policy='omit'" will probably not work.` – kop May 12 '20 at 11:43
  • I also tried to run the 3rd attempt I wrote before, including: `gmodel=Model(model)` `model_values=gmodel.eval(p, t=t)` Now, the model_values is a series full of 1. . This is strange, because when I call the model function assuming that the parameters are known, for example: `relative_fluxNR=model(t=dataOT['time'], Rp=0.000394786, Rs=0.003221125, a=1.86, orb_inclination=89.997*np.pi/180, Rin=0.000394786, Rout=0.000394786, tau=0 )`, the function works properly: [link](https://i.stack.imgur.com/TcmnI.png) – kop May 12 '20 at 11:44
  • Well, the good news is that now you know where to look. Generating NaNs will absolutely kill a fit - you have to find where they are coming from and make sure you do not generate them. Again, I strongly recommend writing your model function so that it can possibly be read. Also, if calling `gmodel.eval(params, t=t)` is different from `gmodel(t, ....)` then your parameters must not have those values. Again, the good news is that you know where to look. – M Newville May 12 '20 at 15:36
  • Based on your recommendation, I redefined function A in order to make the code readable -I'm not using the for loop anymore. I think the problem is that the variable t is not appearing explicitly in the model expression. This expression is the result of the combination of many functions, and it seems impossible to write it in a simpler way (I mean with t appearing in the final expression). Could this be the problem? I've also checked the model values, and there are no Nan values. Note that the model function has the expected behavior in the program, except from the part related with lmfit. – kop May 13 '20 at 20:27
0

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.

kop
  • 21
  • 5