SS2012-StatisticalMethods.git

final version
Parent 40316fc
Michi

Michi commited on 2012-06-06 16:00:53
Zeige 1 geänderte Dateien mit 259 Einfügungen und 0 Löschungen.

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