lmfit 和 scipy curve_fit return 初始猜测作为最佳拟合参数

Question

我想对某些数据拟合一个函数，但我遇到了问题。我尝试使用 lmfit 或 scipy 中的 curve_fit。下面我描述一下问题。

这是我的数据：

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

此外，这里是要拟合数据的模型函数：

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))

其中u1、u2为已知数，待拟合参数为：Rp、Rs、a、orb_inclination、Rin、Rout、tau，它们包含在Agps、Agos、Agis量中。下面是函数 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))

第一次尝试：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)

第二次尝试：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))

所有情况 return 作为初始猜测的最佳拟合参数。我应该怎么做才能使其正常工作？

注意：假设参数已知，模型具有预期的行为（Figure 1），所以我认为模型定义明确，问题与此无关。

如有任何帮助，我们将不胜感激。提前致谢！

Answer 1

我有两个想法，我认为第一个可能是你的罪魁祸首。

在我看来，使用 nan_policy = 'omit' 似乎适用于非常特殊的情况。如果您在尝试拟合时收到错误消息 "nan_policy = 'omit' will probably not work"，那么，它可能无法正常工作。您可以对 NaN 值进行简单检查，以确认该函数针对您的区间输出 NaN 值。
变量的范围很大。尝试提高间隔的最小值。

Answer 2

没有真实的数据和完整的例子，很难猜出哪里出了问题。因此，这将包括一些关于如何解决问题的建议。

首先：由于您正在进行曲线拟合，并且具有模型功能，我建议您从第二个版本开始，使用 lmfit.Model。但是，我建议明确制作一组参数，如：

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())

这本身并不能解决问题，但清晰度很重要。但是，您可以自己调用该模型函数或执行

  gmodel.eval(params, t=t)

并查看它为任何一组参数值实际计算的结果。

其次：您应该谨慎对待跨越多个数量级的拟合问题中的变量。让变量更像阶数 1（或者，在阶数 1.e-6 和 1.e6 之间），然后适当地乘以 1e9 或 1e12 的因数 - 或者只是在值更接近的单位中工作to 1. fitting的数值都是双精度浮点数，参数的相对值很重要

第三：你的模型函数，哎呀。可读性计数。编写一个难以理解的函数对任何人都没有帮助。包括你。我保证你不知道这是做什么的。例如，您可能能够避免循环并仅使用 numpy 的 ufunc-ness，但这是不可能的。需要明确的是，这是不可能分辨的，因为你是这样写的。就像 u1 和 u2 到底应该是什么？真的，这个函数不存在，你写得一团糟，然后出了问题。

所以：写你的模型函数，就好像你希望明年能读到它一样，然后用合理的输入值测试它的计算结果。当它起作用时，合身性也应该起作用。

Answer 3

我通过减少参数数量解决了这个问题。另外，另一个问题是其中一个参数根本不影响拟合。

lmfit 和 scipy curve_fit return 初始猜测作为最佳拟合参数

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

curve-fitting

best-fit-curve

model-fitting

lmfit

scipy-optimize