0

我有如下功能

q = 1 / sqrt( ((1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*(1-o_m)) )
h = 5 * log10( (1+z)*q ) + 43.1601

我有上述方程的实验答案，一旦我必须将一些数据放入上述函数并求解下面的方程

chi=(q_exp-q_theo)**2/err**2  # this function is a sigma, sigma chi from z=0 to z=1.4 (in the data file)

z，err并且q_exp在数据文件（2.txt）中。现在我必须选择一个范围o_m (0.2 to 0.4)并在什么中找到o_m，该chi功能将被最小化。

我的代码是：

from math import *
from scipy.integrate import quad

min = None
l = None
a = None
b = None
c = 0


def ant(z,om,od):
    return 1/sqrt( (1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*o_d )


for o_m in range(20,40,1):
    o_d=1-0.01*o_m
    with open('2.txt') as fp:
        for line in fp:
            n = list( map(float, line.split()) )
            q = quad(ant,n[0],n[1],args=(o_m,o_d))[0]
            h = 5.0 * log10( (1+n[1])*q ) + 43.1601
            chi = (n[2]-h)**2 / n[3]**2
            c = c + chi
        if min is None or min>c:
            min = c
            l = o_m                            
        print('chi=',q,'o_m=',0.01*l)

n[1], n[2], n[3],n[4]是z1, z2,q_exp和err, 分别在数据文件中。和是积分范围z1。z2我需要您的帮助，感谢您的时间和关注。请不要评价负值。我需要你的答案。

4

2 回答 2

1

这是我对问题的理解。首先我通过以下代码生成一些数据

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协方差矩阵是简化形式，需要重新调整。

于 2017-08-04T14:35:06.240 回答
0

不自觉地，这个问题与我的另一个问题重叠。正确答案是：

 from math import *
 import numpy as np
 from scipy.integrate import quad
 min=l=a=b=chi=None
 c=0
 z,mo,err=np.genfromtxt('Union2.1_z_dm_err.txt',unpack=True)
 def ant(z,o_m):            #0.01*o_m  is steps of o_m
     return 1/sqrt(((1+z)**2*(1+0.01*o_m*z)-z*(2+z)*(1-0.01*o_m)))
 for o_m in range(20,40):
     c=0
     for i in range(len(z)):
         q=quad(ant,0,z[i],args=(o_m,))[0]     #Integration o to z
         h=5*log10((1+z[i])*(299000/70)*q)+25     #function of dL
         chi=(mo[i]-h)**2/err[i]**2               #chi^2 test function
        c=c+chi
        l=o_m
        print('chi^2=',c,'Om=',0.01*l,'OD=',1-0.01*l)
于 2017-08-14T06:58:09.840 回答