c64b192097f74ba5493b0faadd50feaf2c105d25
Michi add sheet 6

Michi authored 11 years ago

1) #Sheet 6
2) 
3) #Problem1
4) 
5) #setwd("C:/Users/Studium/Desktop/Statistische Methoden/Sheet 6")
6) input=scan(file = "sn_data_riess.dat", what = list(character(), double(), double(), double()), skip=1, multi.line=FALSE)
7) data=cbind(input[[2]], input[[3]], input[[4]])
8) 
9) #Hubble expansion in Dark Energy models
10) #we assume flat Universe
11) H0=72.0 #km/s/Mpc
12) c=3*10^5 #km/s
13) N=186 #number of data
14) 
15) #################Model_A########################
16) #Hubble law for LCDM
17) Hubble_A=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))
18) 
19) #H^(-1)
20) Hinv_A=function(zint,Omegam) 1.0/Hubble_A(zint,Omegam)
21) 
22) #Luminosity distance
23) dL_A=function(z,Omegami){ 
24) 	dLsol=z
25) 	for (i in (1:length(z))){
26) 		zarg=z[i]
27) 		I=integrate(Hinv_A,0.0,zarg,Omegam=Omegami)
28) 		dLsol[i] = c*(1+zarg)*I$value
29) 		}
30) 	return(dLsol)
31) 	}
32) 
33) #Magnitude
34) mag_A=function(z,Omegam,M) M+5.0*log10(H0*dL_A(z,Omegam))
35) 
36) # Define function chi2 
37) chi2_A=function(Omegam,M){
38) chi2sol=0.0
39) for (i in (1:N))  
40)                {          
41) 		chi2sol=chi2sol+((data[i,2]-mag_A(data[i,1],Omegam,M))^2)/(data[i,3])^2
42) 	       }
43) 	return(chi2sol)
44) }
45) 
46) ###############Model_B####################
47) #Hubble law for wCDM
48) Hubble_B=function(z,Omegam,w) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam)*(1+z)^(3*(1+w)))
49) 
50) #H^(-1)
51) Hinv_B=function(zint,Omegam, w) 1.0/Hubble_B(zint,Omegam, w)
52) 
53) 
54) #Luminosity distance
55) dL_B=function(z,Omegami, wi){ 
56) 	dLsol=z
57) 	for (i in (1:length(z))){
58) 		zarg=z[i]
59) 		I=integrate(Hinv_B,0.0,zarg,Omegam=Omegami, w=wi)
60) 		dLsol[i] = c*(1+zarg)*I$value
61) 		}
62) 	return(dLsol)
63) 	}
64) 
65) #Magnitude
66) mag_B=function(z,Omegam,M, w) M+5.0*log10(H0*dL_B(z,Omegam, w))
67) 
68) # Define function chi2 
69) chi2_B=function(Omegam,M ,w){
70) chi2sol=0.0
71) for (i in (1:N))  
72)                {          
73) 		chi2sol=chi2sol+((data[i,2]-mag_B(data[i,1],Omegam,M,w))^2)/(data[i,3])^2
74) 	       }
75) 	return(chi2sol)
76) }
77) 
78) 
79) #i) Omegam=0,3, w=m{-2.0, -1.5, -1.0, -0.5}
80) #plot(data[,1], dL_B(data[,1], 0.3, -2.0),type="p", col="black", pch=19, cex=.6)
81) #points(data[,1], dL_B(data[,1], 0.3, -1.5),type="p", col="blue", pch=19, cex=.6)
82) #points(data[,1], dL_B(data[,1], 0.3, -1.0),type="p", col="green", pch=19, cex=.6)
83) #points(data[,1], dL_B(data[,1], 0.3, -0.5),type="p", col="red", pch=19, cex=.6)
84) 
85) #windows()
86) #ii) w=-1, Omegam={0.1, 0.2, 0.3, 0.4, 0.5}
87) #plot(data[,1], dL_B(data[,1], 0.1, -1.0),type="p", col="black", pch=19, cex=.6)
88) #points(data[,1], dL_B(data[,1], 0.2, -1.0),type="p", col="blue", pch=19, cex=.6)
89) #points(data[,1], dL_B(data[,1], 0.3, -1.0),type="p", col="green", pch=19, cex=.6)
90) #points(data[,1], dL_B(data[,1], 0.4, -1.0),type="p", col="red", pch=19, cex=.6)
91) #points(data[,1], dL_B(data[,1], 0.5, -1.0),type="p", col="yellow", pch=19, cex=.6)
92) 
93) 
94) #Problem 2
95) 
96) 
97) #windows()
98) #i) Model_A: Omegam variiert von 0.1 bis 0.5
99) #plot(data[,1], mag_A(data[,1], 0.1, 16),type="p", col="black", pch=19, cex=.6)
100) #points(data[,1], mag_A(data[,1], 0.2, 16),type="p", col="blue", pch=19, cex=.6)
101) #points(data[,1], mag_A(data[,1], 0.3, 16),type="p", col="green", pch=19, cex=.6)
102) #points(data[,1], mag_A(data[,1], 0.4, 16),type="p", col="red", pch=19, cex=.6)
103) #points(data[,1], mag_A(data[,1], 0.5, 16),type="p", col="yellow", pch=19, cex=.6)
104) 
105) #ii)
106) #windows()
107) # Model_B: Omegam variiert von 0.1 bis 0.5 bei w=-1.0
108) #plot(data[,1], mag_B(data[,1], 0.1, 16, -1.0),type="p", col="black", pch=19, cex=.6)
109) #points(data[,1], mag_B(data[,1], 0.2, 16, -1.0),type="p", col="blue", pch=19, cex=.6)
110) #points(data[,1], mag_B(data[,1], 0.3, 16, -1.0),type="p", col="green", pch=19, cex=.6)
111) #points(data[,1], mag_B(data[,1], 0.4, 16, -1.0),type="p", col="red", pch=19, cex=.6)
112) #points(data[,1], mag_B(data[,1], 0.5, 16, -1.0),type="p", col="yellow", pch=19, cex=.6)
113) 
114) #windows()
115) # Model_B: w variiert von -2.0 bis -0.5 bei Omegam=0.3
116) #plot(data[,1], mag_B(data[,1], 0.3, 16, -2.0),type="p", col="black", pch=19, cex=.6)
117) #points(data[,1], mag_B(data[,1], 0.3, 16, -1.5),type="p", col="blue", pch=19, cex=.6)
118) #points(data[,1], mag_B(data[,1], 0.3, 16, -1.0),type="p", col="green", pch=19, cex=.6)
119) #points(data[,1], mag_B(data[,1], 0.3, 16, -0.5),type="p", col="red", pch=19, cex=.6)
120) 
121) #Problem 3
122) #Model_A: chi2 for Omegam={0.0,1.0} and M={15.5, 16.5}
123) values_A=10
124) Omegam_vec_A=seq(0.0, 1.0, by = 1/values_A)
125) M_vec_A=seq(15.5, 16.5, by=1/values_A)
126) twodim=array(0, dim=c(values_A+1,values_A+1))
127) for (i in (0:values_A+1)){
128) 	for (j in (0:values_A+1)){
129) 		twodim[i,j]=chi2_A(Omegam_vec_A[i], M_vec_A[j])
130) 	}
131) }
132) 
133) #Model_B: chi2 for Omegam={0.0,1.0} and M={15.5, 16.5}, w={-2.0,-0.5}
134) values_B=5
135) Omegam_vec_B=seq(0.0, 1.0, by = 1/values_B)
136) M_vec_B=seq(15.5, 16.5, by=1/values_B)
137) w_vec=seq(-2.0, -0.5, by=1/(values_B/1.5))
138) 
139) threedim=array(0, dim=c(values_B+1,values_B+1, values_B+1))
140) for (i in (0:values_B+1)){
141) 	for (j in (0:values_B+1)){
142) 		for (k in (0:values_B+1)){
143) 			threedim[i,j,k]=chi2_B(Omegam_vec_B[i], M_vec_B[j], w_vec[k])
144) 		}
145) 	}
146) }
147) #Find best fit parameters: where is the minimum of chi2?
148) #Model_A: find minimum
149) min_A=twodim[1,1]
150) min_A_pos=array(0, dim=c(2))
151) for (i in (0:values_A+1)){
152) 	for (j in (0:values_A+1)){
153) 		if (twodim[i,j] < min_A){
154) 			min_A=twodim[i,j]
155) 			min_A_pos[1]=i
156) 			min_A_pos[2]=j
157) 		}
158) 	}
159) }
160) twodim
161) out1=c("chi2min",min_A)
162) out2=c("Omegam_min",Omegam_vec_A[min_A_pos[1]])
163) out3=c("M_min",M_vec_A[min_A_pos[2]])                                
164) print(out1)
165) print(out2)
166) print(out3)
167) 
168) #Model_B: find minimum
169) min_B=threedim[1,1,1]
170) min_B_pos=array(0, dim=c(3))
171) for (i in (0:values_B+1)){
172) 	for (j in (0:values_B+1)){
173) 		for (k in (0:values_B+1)){
174) 		
175) 			if (threedim[i,j,k] < min_B){
176) 				min_B=threedim[i,j,k]
177) 				min_B_pos[1]=i
178) 				min_B_pos[2]=j
179) 				min_B_pos[3]=k
180) 			}
181) 		}
182) 	}
183) }
184) threedim
185) out1=c("chi2min",min_B)
186) out2=c("Omegam_min",Omegam_vec_B[min_B_pos[1]])
187) out3=c("M_min",M_vec_B[min_B_pos[2]]) 
188) out4=c("w_min",w_vec[min_B_pos[3]])                               
189) print(out1)
190) print(out2)
191) print(out3)
192) print(out4)
193) 
194) #Problem4
195) 
196) #Model_A posterior function  
197) posterior_A=array(0, c(values_A+1, values_A+1)) 
198) for (i in (0:values_A)){
199) 	for (j in (0:values_A)){
200) 		posterior_A[i,j]=exp(-0.5*(twodim[i,j]-min_A))
201) 	}
202) }                    
203)   
204) #Model_B posterior function                        
205) posterior_B=array(0, c(values_B+1, values_B+1, values_B+1))
206) for (i in (0:values_B)){
207) 	for (j in (0:values_B)){
208) 		for (k in (0:values_B)){
209) 			posterior_B[i,j,k]=exp(-0.5*(threedim[i,j,k]-min_B))
210) 		}
211) 	}
212) }    
213) 
214) give_function=function(i,j) f=splinefun(w_vec, posterior_B[i,j,])
215) Myintegrate = function(F,LOW,UP){
216)     pI = integrate(F,LOW,UP)
217)     return(pI$value)
218) }
219) 
220) integration_B=array(0, c(values_B+1, values_B+1))
221) for (i in (1:values_B)){
222) 	for (j in (1:values_B)){
223) 		I=Myintegrate(give_function(i,j),-2.0,-0.5)
224) 		integration_B[i,j]=I
225) 	}
226) }
227) 
228) 
229) contour(Omegam_vec_A,M_vec_A,posterior_A,drawlabels=FALSE,xlab='Omegam',ylab='M',xlim=c(0,1),ylim=c(15.5,16.5),col = "red")
230) contour(Omegam_vec_B,M_vec_B,integration_B,add = TRUE)
231) 
232) # Marginalize over M
Michi add 6.5

Michi authored 11 years ago

233) 
234) 
235) give_function_MA=function(i) f=splinefun(M_vec_A, posterior_A[i,])
236) give_function_MB=function(i,ii) f=splinefun(M_vec_B, posterior_B[i,,ii])
237) 
238) 
239) integration_MA=array(0, c(values_A+1))
240) for (i in (1:values_A)){
241) 		I=Myintegrate(give_function_MA(i),15.5,16.5)
242) 		integration_MA[i]=I
243) 	}
244) 
Michi dimension mismatch in 6.5

Michi authored 11 years ago

245) integration_MB=array(0, c(values_B, values_B))
Michi add 6.5

Michi authored 11 years ago

246) for (i in (1:values_B)){
247) 	for (ii in (1:values_B)){
248) 		I=Myintegrate(give_function_MB(i,ii),15.5,16.5)
249) 		integration_MB[i,ii]=I
250) 	}
251) }
252) 
Michi dimension mismatch in 6.5

Michi authored 11 years ago

253) contour(Omegam_vec_B,w_vec,integration_MB)