Тест хи-квадрат для ограничения параметра

У меня важный вопрос об использовании критерия chi ^ 2 для ограничения параметра в космологии. Я ценю вашу помощь. Пожалуйста, не ставьте этому вопросу отрицательную оценку (этот вопрос важен для меня).

Предположим, у нас есть файл данных (data.txt), содержащий 600 данных, и этот файл данных имеет 3 столбца, первый столбец - красное смещение (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

Этот вопрос не повторяется и очень важен для каждого человека, использующего ци ^ 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)

Я думаю, что мой код правильный, но он не дает мне правильного ответа. Спасибо, и я ценю ваше время и внимание.


python-3.x python python-2.7 data-fitting chi-squared
person Ethan    schedule 11.08.2017    source источник
comment
Этот вопрос не по теме в Stackoverflow, поскольку он не касается конкретной проблемы программирования. Однако вы можете взглянуть на scipy.optimize для численная оптимизация. Если вы хотите решить свою проблему аналитически, вам следует рассмотреть возможность использования метода максимального правдоподобия.   -  person MaxPowers    schedule 11.08.2017
comment
@maxpowers, это как раз проблема, я написал код, но он не дает мне правильного ответа. Спасибо за помощь   -  person Ethan    schedule 11.08.2017
comment
Прикрепил код, может поможет. Спасибо. Пожалуйста, не ставьте мне отрицательную оценку. Это несправедливо, друг мой.   -  person Ethan    schedule 11.08.2017

Ответы (1)


arrow_upward
2
arrow_downward

Это правильный ответ:

 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)
person Community    schedule 14.08.2017