chi-square testing for constraining a parameter

2019-06-05 18:11发布

问题:

I have an important question about the use of chi^2 test to constrain a parameter in cosmology. I appreciate your help. Please do not give this question negative rate (this question is important to me).

Assume we have a data file(data.txt) concluding 600 data and this data file has 3 columns, first column is redshift(z), second column is observational dL(m_obs) and third column is error(err). As we know chi^2 function is

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

All of thing that we must calculate is putting "z" from given data file into our function in "m_theo" for all 600 data and calculate chi^2. Now in "m_thoe" we have a free parameter(o_m) and we have to find its value in which the chi^2 reach its minimum value.

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

This question is not repetitive and is very important to every body using chi^2 specially for cosmologists and physicists. How to find minimized chi^2 and relative 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)

I think my code is correct but it doesn't give me a proper answer Thank you and I appreciate your time and your attention.

回答1:

This is the correct answer:

 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)