2017-08-11 43 views
1

我有一个关于使用chi^2测试来限制宇宙学参数的重要问题。我感谢您的帮助。请不要给这个问题负面的利率(这个问题对我很重要)。假设我们有一个结束600个数据的数据文件(data.txt),这个数据文件有3列,第一列是红移(z),第二列是观测dL(m_obs),第三列是错误(err) 。正如我们所知道志^ 2的功能是约束参数的卡方检验

chi^2=(m_obs-m_theo)**2/err**2 #chi^2=sigma((m_obs-m_theo)**2/err**2) from 1 to N=600 

的,我们必须计算从给定的数据文件中把“Z”到我们的功能“m_theo”的所有600个数据,并计算卡^ 2的事情全部。现在在“m_thoe”中我们有一个自由参数(o_m),我们必须找到其中chi^2达到其最小值的值。

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

这个问题不重复使用是每一个身体非常重要的阴气^ 2专为宇宙学家和物理学家。 如何找到最小化chi^2和相对o_m?

from math import * 
import numpy as np 
from scipy.integrate import quad 
min=l=a=b=chi=None 
c=0 #for Sigma or summation chi^2 terms in c=c+chi for first term 
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(24,35,1): #arbitrary range of o_m 
############## opening data file containing 580 dataset 
    with open('data.txt') as f: 
     for i, line in enumerate(f): # 
      n= list(map(float, line.split())) # 
      for i in range(1): 
##############    
       q=quad(ant,0,n[1],args=(o_m,))[0] #Integration o to z, z=n[1] 
       h=5*log10((1+n[1])*(299/70)*q)+25 #function of dL 
       chi=(n[2]-h)**2/n[3]**2  #chi^2 test function 
       c=c+chi   #sigma from 1 to N of chi^2 and N=580 
     if min is None or min>c: 
      min=c 
      print(c,o_m) 

我觉得我的代码是正确的,但它并没有给我一个正确的答案 谢谢你,我很欣赏你的时间和注意力。

+0

这个问题是题外话#2,因为它是不是一个独特的规划问题。但是,您可能希望查看[scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html)进行数值优化。如果你想解析你的问题,你应该考虑使用[最大似然法](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation)。 – MaxPowers

+0

@maxpowers这完全是一个问题,我写了一个代码,但它没有给我正确的答案。谢谢你的帮助 – Ethan

+1

我附上了代码,可能有帮助。谢谢。请不要给我负面的利率。我的朋友不公平。 – Ethan

回答

1

这是正确的答案:

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)