scipy的Odeint和Curve fit会发出警告并给出错误的结果

如何解决scipy的Odeint和Curve fit会发出警告并给出错误的结果

我对odeint和曲线拟合的警告感到困扰。所以我想做的是：

我的第一个问题是，对于3个数据集，曲线拟合和odeint反复发出如下警告（总共6条警告），但是curve_fit确实给出了看似正确的结果。

828：OptimizeWarning：无法估计参数的协方差

247：ODEintWarning：在此调用上完成了过多的工作（也许是错误的Dfun类型）。以full_output = 1运行以获取定量信息。

**第二个问题是积分曲线，通过多次执行EXACT SAME代码，它会产生不同的结果。似乎有一段重复，在执行了几次之后，它给出了正确的曲线，但是在下一次执行时，它再次变得不正确，依此类推。也许是不稳定的问题？**

import math
import numpy as np
import pathlib as pl
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.integrate import odeint

def loadData(path):    
    with open(path,"r") as fid:
        res=np.loadtxt(fid,comments="#")
    return res

def modeleeq(a,N,p,we,wp):
    res=(we*a/p[0])**p[1] + (wP*a/p[2])**p[3]
    return res

#get path
chemin=pl.Path(input("Paste the path to the directory containing data files\n"))

#get we and wp in array 3X2
wewp=loadData(chemin/'wewp.dat').transpose()

#plotting
plt.style.use('seaborn')
fig,(ax1,ax2) = plt.subplots(1,2)
colors = ['#79ccff','#f78db4','#a07ffb']

files = pl.Path(chemin).glob("a_dadN*")
para=[]
for i,f in enumerate(files):
    res=loadData(f)
    a=res[:,0]
    dadN=res[:,1]
    we=wewp[i,0]
    wp=wewp[i,1]
    
    #exp data    
    ax1.scatter(a,dadN,c=colors[i],marker="x",label = f"test {i+1}")
    
    #evaluation of parameters
    modele=lambda a,*p: (we*a/p[0])**p[1] + (wP*a/p[2])**p[3]  #p=[ge,me,gp,mp]
    p,pcov= curve_fit(modele,a,p0 = [2e5,1e5,0])
    para.append(p)
    ax1.plot(a,modele(a,*p),label=f"test {i+1} identification")
    
    #intergration
    a0=a.min() #initial condition
    N = np.linspace(0,4000)
    aitgr=odeint(modeleeq,a0,args=(p,wp))
    ax2.plot(N,aitgr,label = f"test {i+1} integration")
#the following code is just for adding titles and print the identified paras
#so I won't put them here

非常感谢大家！

解决方法

我认为这是一些问题。首先，很明显，这两个加数非常相似。因此，数据确实无法区分两者，这确实可能发生。扩展是第二个问题。具有跨越数量级的拟合参数通常不是一个好主意。实际上，伽马只是重新缩放W。因此，人们可以轻松地将函数重写为( f1 * a )**e1，因为它适合f1并通过知道gamma来计算f1=W/gamma。另一个问题是出现负数的负幂的可能性。因此，应该使用abs( f1 )或f1**2。考虑到这一点，我修改了代码并得到以下结果。从第二和第三种情况的拟合结果中，可以看出f1=0或指数几乎相等。在这种情况下，通常无法确定协方差矩阵。

最后，在积分/绘制a(N)时，微分方程的类型为da / dN = a**k，因此我们有da / a**k = dN积分将得到a**(-k+1) = N+N0，所以a（N） = 1 /（N + N0）**（1 /（k-1））。由于初始斜率为N0 < 0为正，并且函数在N0处发散。超出此值的数值积分没有意义。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.integrate import odeint

def loadData(path):    
    with open(path,"r") as fid:
        res=np.loadtxt(fid,comments="#")
    return res


def simplified_model(a,N,f1,e1,f2,e2 ):
    dadN = ( np.fabs( f1 )  * a )**e1  + ( np.fabs( f2 ) * a )**e2
    return dadN

plt.style.use('seaborn')
fig = plt.figure()
ax1 = fig.add_subplot( 1,2,1)
ax2 = fig.add_subplot( 1,2)

colors = ['#79ccff','#f78db4','#a07ffb']

para = list()
for i in range( 3 ):
    res=loadData( "a_dadN_{}.dat".format( i + 1 ) )
    a = res[ ::,0 ]
    dadN = res[ ::,1 ]
    # ~we = wewp[i,0 ]
    # ~wp = wewp[ i,1 ]

    def model_wrapper(a,ge,me,gp,mp ):
        return simplified_model(a,mp )

    #exp data    
    ax1.scatter( a,dadN,c=colors[i],marker="x",label = f"test {i+1}" )

    print(i)

    al = np.linspace( min(a),max(a),50 )
    guess = [ .2 + i,2.0,2.2 ]
    dal = np.fromiter( ( model_wrapper( av,*guess ) for av in al ),np.float )
    ax1.plot( al,dal,color=colors[i],ls=':')
    p,pcov= curve_fit(model_wrapper,a,p0 = guess,maxfev=100000 )
    print( p )
    dalf = np.fromiter( ( model_wrapper( av,*p ) for av in al ),dalf,color=colors[i])

    #intergration
    a0 = a.min() #initial condition
    N = np.linspace(0,4000,100000)
    aitgr=odeint( simplified_model,a0,args=tuple(p) )
    ax2.plot( N,aitgr,label = f"test {i+1} integration")
    ax2.set_ylim([1e-4,1e-2])
    ax2.set_yscale("log")

plt.show()

提供