问题描述
我对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()
提供