SS2012-StatisticalMethods.git

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#Sheet 6
#Problem1
#setwd("C:/Users/Studium/Desktop/Statistische Methoden/Sheet 6")
input=scan(file = "sn_data_riess.dat", what = list(character(), double(), double(), double()), skip=1, multi.line=FALSE)
data=cbind(input[[2]], input[[3]], input[[4]])
#Hubble expansion in Dark Energy models
#we assume flat Universe
H0=72.0 #km/s/Mpc
c=3*10^5 #km/s
N=186 #number of data
#################Model_A########################
#Hubble law for LCDM
Hubble_A=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))
#H^(-1)
Hinv_A=function(zint,Omegam) 1.0/Hubble_A(zint,Omegam)
#Luminosity distance
dL_A=function(z,Omegami){ 
    dLsol=z
    for (i in (1:length(z))){
        zarg=z[i]
        I=integrate(Hinv_A,0.0,zarg,Omegam=Omegami)
        dLsol[i] = c*(1+zarg)*I$value
        }
    return(dLsol)
    }
#Magnitude
mag_A=function(z,Omegam,M) M+5.0*log10(H0*dL_A(z,Omegam))
# Define function chi2 
chi2_A=function(Omegam,M){
chi2sol=0.0
for (i in (1:N))  
               {          
        chi2sol=chi2sol+((data[i,2]-mag_A(data[i,1],Omegam,M))^2)/(data[i,3])^2
           }
    return(chi2sol)
}
###############Model_B####################
#Hubble law for wCDM
Hubble_B=function(z,Omegam,w) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam)*(1+z)^(3*(1+w)))
#H^(-1)
Hinv_B=function(zint,Omegam, w) 1.0/Hubble_B(zint,Omegam, w)
 
 
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX