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

问题描述 投票：0回答：2

我想对某些数据使用函数，但遇到了问题。我尝试从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):  
    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

第一次尝试：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 Model
gmodel=Model(model)
result = gmodel.fit(y, t=t, 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)
print(result.fit_report())

第三次尝试：lmfit最小化器

from lmfit import Parameters, Minimizer, report_fit
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=4.5*10**9, min= 0,max=8*10**9)
p.add('Rs' ,  value=4.3*10**10, min= 0,max=9*10**10)
p.add('a',   value=1.4*10**13, min= 0,max= 4*10**13)
p.add('orb_inclination',  value=1, min= 0,max=4)
p.add('Rin',  value=4.5*10**9, min= 0,max=8*10**9)
p.add('Rout',  value=4.5*10**9, min= 0,max=23*10**9)
p.add('tau',  value=0,min= 0,max=2)

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

# first solve with Nelder-Mead
out1 = mini.minimize(method='Nelder')

# then solve with Levenberg-Marquardt using the
# Nelder-Mead solution as a starting point
out2 = mini.minimize(method='leastsq', params=out1.params)
print(report_fit(out2))

所有情况均以最合适的参数返回初始猜测。我应该怎么做才能使其正常工作？

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

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

curve-fitting scipy-optimize model-fitting lmfit best-fit-curve
2个回答
0
投票

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

  1. 在我看来，使用nan_policy ='省略'可以在非常特殊的情况下工作。如果在尝试拟合时收到错误消息“ nan_policy ='omit'可能不起作用”，那么可能就行不通了。您可以对NaN值进行简单检查，以确认该函数为您的间隔输出NaN值。
  2. 变量的边界为huge。尝试提高间隔的最小值。
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到底应该是什么？确实，此功能不存在，您编写了一个完整的烂摊子，然后出了点问题。

因此：编写模型函数，就像您希望明年阅读它一样，然后使用合理的输入值测试它的计算结果。当这种方法起作用时，拟合也应该起作用。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.