1

我想为一些数据拟合一个函数，但我遇到了一个问题。我尝试使用 scipy 中的 lmfit 或 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))

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

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

任何帮助，将不胜感激。先感谢您！

4

3 回答 3

0

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

  1. 在我看来，使用 nan_policy = 'omit' 似乎适用于非常具体的情况。如果您在尝试适应时收到错误消息“nan_policy = 'omit' 可能不起作用”，那么它可能不起作用。您可以对 NaN 值进行简单检查，以确认该函数为您的时间间隔输出 NaN 值。
  2. 变量的界限是巨大的。尝试提高间隔的最小值。
于 2020-05-11T15:06:40.800 回答
0

如果没有真实数据和完整示例，很难猜测可能出了什么问题。因此，这将包括一些关于如何解决问题的建议。

第一：因为你正在做曲线拟合，并且有一个模型功能，我建议从你的第二个版本开始，使用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 的因子 - 或者仅以值接近 1 的单位工作。数字拟合的都是双精度浮点数，参数的相对值很重要。

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

所以：编写你的模型函数，就像你希望明年阅读它一样，然后用合理的输入值测试它的计算结果。当这有效时，合身也应该有效。

于 2020-05-12T02:45:49.153 回答
0

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

于 2020-05-19T05:17:23.360 回答