python - 用于约束参数的卡方检验

标签 python python-2.7 python-3.x data-fitting chi-squared

我有一个关于使用 chi^2 检验来约束宇宙学参数的重要问题。我感谢您的帮助。请不要给这个问题负分(这个问题对我来说很重要)。

假设我们有一个包含 600 个数据的数据文件 (data.txt),该数据文件有 3 列,第一列是 redshift(z),第二列是观测 dL(m_obs),第三列是列有错误(err)。我们知道 chi^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 个数据并计算 chi^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

这个问题不是重复的,对于每个使用 chi^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)

我认为我的代码是正确的,但它没有给我正确的答案。 谢谢您,感谢您的时间和关注。

最佳答案

这是正确答案:

 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)

关于python - 用于约束参数的卡方检验,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/45631732/

相关文章:

python datetime.utcnow 不显示正确的时间戳

python - 如何使用 subprocess.Popen 更新目录以在 Python 中查找新文件

python - django-import-export 覆盖不起作用

python - 枕头 : File not found

python - openerp错误AttributeError : 'int' object has no attribute 'iteritems'

Python:用单词列表替换句子中的一个单词,并将新句子放在 pandas 的另一列中

python-3.x - 如何通过boto库为EMR集群选项配置 "Use AWS Glue Data Catalog for table metadata"?

python - 如何解析 python (spyder) 中的参数?

python - 亚马逊网络服务 (AWS) 上的 Django

Python - 从 lxml xpath 获取类