这里是我的问题的理解。 首先,我通过下面的代码
import numpy as np
from scipy.integrate import quad
from random import random
def boxmuller(x0,sigma):
u1=random()
u2=random()
ll=np.sqrt(-2*np.log(u1))
z0=ll*np.cos(2*np.pi*u2)
z1=ll*np.cos(2*np.pi*u2)
return sigma*z0+x0, sigma*z1+x0
def q_func(z, oM, oD):
den= np.sqrt((1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD))
return 1.0/den
def h_func(z,q):
out = 5 * np.log10((1.0 + z) * q) + .25#43.1601
return out
def q_Int(z1,z2,oM,oD):
out=quad(q_func, z1,z2,args=(oM,oD))
return out
ooMM=0.3
ooDD=1.0-ooMM
dataList=[]
for z in np.linspace(.3,20,60):
z1=.1+.1*z*.01*z**2
z2=z1+3.0+.08+z**2
q=q_Int(z1,z2,ooMM,ooDD)[0]
h=h_func(z,q)
sigma=np.fabs(.01*h)
h=boxmuller(h,sigma)[0]
dataList+=[[z,z1,z2,h,sigma]]
dataList=np.array(dataList)
np.savetxt("data.txt",dataList)
我将后装配以如下方式
import matplotlib
matplotlib.use('Qt5Agg')
from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.optimize import leastsq
def q_func(z, oM, oD):
den= np.sqrt((1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD))
return 1.0/den
def h_func(z,q):
out = 5 * np.log10((1.0 + z) * q) + .25#43.1601
return out
def q_Int(z1,z2,oM,oD):
out=quad(q_func, z1,z2,args=(oM,oD))
return out
def residuals(parameters,data):
om,od=parameters
zList=data[:,0]
yList=data[:,3]
errList=data[:,4]
qList=np.fromiter((q_Int(z1,z2, om,od)[0] for z1,z2 in data[ :,[1,2] ]), np.float)
hList=np.fromiter((h_func(z,q) for z,q in zip(zList,qList)), np.float)
diffList=np.fromiter(((y-h)/e for y,h,e in zip(yList,hList,errList)), np.float)
return diffList
dataList=np.loadtxt("data.txt")
###fitting
startGuess=[.4,.8]
bestFitValues, cov,info,mesg, ier = leastsq(residuals, startGuess , args=(dataList,),full_output=1)
print bestFitValues,cov
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(dataList[:,0],dataList[:,3],marker='x')
###fitresult
fqList=[q_Int(z1,z2,bestFitValues[0], bestFitValues[1])[0] for z1,z2 in zip(dataList[:,1],dataList[:,2])]
fhList=[h_func(z,q) for z,q in zip(dataList[:,0],fqList)]
ax.plot(dataList[:,0],fhList,marker='+')
plt.show()
给输出
>>[ 0.31703574 0.69572673]
>>[[ 1.38135263e-03 -2.06088258e-04]
>> [ -2.06088258e-04 7.33485166e-05]]
和图形 注生成一些数据,对于leastsq
协方差矩阵是简化形式,需要重新调整比例。
1:什么是你的问题? 2:请分享一些最小数据集。 3:为什么'ant()'有'o_m'和'o_d'而上面的'q'只有'o_m'。 – mikuszefski
为什么不使用'scipy.optimize.leastsq'?顺便说一句,如果数据不是太大,只需在开始时加载一次,也许用'numpy.loadtxt()' – mikuszefski
此外,还有一些错别字;在最后的'print()'中用'''打开但用'''关闭,你可能想打印'c'而不是'q',你应该命名,例如'c2',因为它实际上是方形的。这样可以避免混淆,有些缩进看起来是错误的。对输入pythonic的评论:'min == None'可以工作,但'min is None'看起来更好,甚至可能是'if not min',你真的想和'h'为'C'? – mikuszefski