SS2012-StatisticalMethods.git

add sheet 6
Parent 8505c07
Michi

Michi commited on 2012-06-03 16:39:58
Zeige 2 geänderte Dateien mit 423 Einfügungen und 0 Löschungen.

6/sn_data_riess.dat
Zeige Datei @ 052b6b9
... ... 
@@ -0,0 +1,187 @@
1 
+#sample  z	   moduli   sigma
2 
+gold     0.0400    36.38    0.19
3 
+gold     0.050     36.84    0.21
4 
+gold     0.0307    35.90    0.20
5 
+gold     0.0560    37.31    0.18
6 
+gold     0.0331    35.54    0.20
7 
+gold     0.0141    34.13    0.25
8 
+gold     0.0460    36.35    0.21
9 
+gold     0.0265    35.64    0.20
10 
+gold     0.101     38.73    0.20
11 
+gold     0.075     37.77    0.19
12 
+gold     0.061     37.30    0.22
13 
+gold     0.0141    34.12    0.25
14 
+gold     0.0262    35.06    0.24
15 
+gold     0.0430    36.53    0.19
16 
+gold     0.0450    36.97    0.18
17 
+gold     0.036     36.17    0.19
18 
+gold     0.058     37.13    0.19
19 
+gold     0.063     37.67    0.19
20 
+gold     0.0186    34.96    0.22
21 
+gold     0.079     37.94    0.18
22 
+gold     0.088     38.07    0.28
23 
+gold     0.0178    34.70    0.23
24 
+gold     0.071     37.78    0.19
25 
+gold     0.0251    35.09    0.21
26 
+gold     0.052     37.16    0.18
27 
+gold     0.0286    35.53    0.21
28 
+gold     0.0490    36.90    0.20
29 
+gold     0.050     37.08    0.19
30 
+gold     0.0180    34.29    0.23
31 
+silver   0.089     38.50    0.17
32 
+silver   0.051     36.67    0.16
33 
+gold     0.0244    35.09    0.20
34 
+gold     0.0290    35.70    0.19
35 
+gold     0.0161    34.50    0.24
36 
+gold     0.0360    36.01    0.20
37 
+silver   0.0116    32.96    0.29
38 
+gold     0.478     42.48    0.23
39 
+silver   0.053     37.17    0.15
40 
+silver   0.230     40.44    0.46
41 
+silver   0.300     40.76    0.60
42 
+silver   0.067     37.54    0.34
43 
+gold     0.450     42.13    0.21
44 
+gold     0.480     42.37    0.20
45 
+gold     0.615     42.85    0.23
46 
+gold     0.400     42.04    0.19
47 
+silver   0.498     43.21    0.24
48 
+silver   0.465     42.81    0.22
49 
+gold     0.0490    36.52    0.20
50 
+gold     0.0219    34.70    0.22
51 
+gold     0.3880    42.07    0.19
52 
+gold     0.0152    34.11    0.25
53 
+gold     0.0276    35.90    0.20
54 
+gold     0.425     41.70    0.40
55 
+gold     0.620     43.11    0.30
56 
+gold     0.570     42.81    0.25
57 
+gold     0.300     41.01    0.25
58 
+gold     0.380     42.02    0.22
59 
+silver   0.160     39.08    0.40
60 
+silver   0.240     40.68    0.43
61 
+gold     0.430     42.33    0.34
62 
+silver   0.0247    35.33    0.25
63 
+gold     0.124     39.20    0.22
64 
+gold     0.0165    33.82    0.27
65 
+gold     0.0167    34.21    0.23
66 
+gold     0.0348    36.17    0.19
67 
+silver   0.490     42.58    0.19
68 
+silver   0.450     42.58    0.19
69 
+gold     0.828     43.96    0.46
70 
+gold     0.495     42.25    0.19
71 
+silver   0.570     42.77    0.19
72 
+gold     0.0132    34.02    0.26
73 
+gold     0.580     43.04    0.21
74 
+gold     0.526     42.56    0.18
75 
+gold     0.172     39.79    0.18
76 
+gold     0.180     39.98    0.18
77 
+gold     0.472     42.46    0.19
78 
+gold     0.430     41.99    0.18
79 
+gold     0.657     43.27    0.20
80 
+gold     0.0166    34.54    0.23
81 
+gold     0.450     42.10    0.23
82 
+gold     0.320     41.45    0.18
83 
+gold     0.581     42.63    0.19
84 
+gold     0.440     42.57    0.40
85 
+gold     0.508     41.64    0.35
86 
+gold     0.416     42.10    0.19
87 
+gold     0.830     43.85    0.19
88 
+gold     0.579     42.86    0.19
89 
+gold     0.420     41.76    0.23
90 
+gold     0.518     42.83    0.30
91 
+gold     0.334     40.92    0.30
92 
+silver   0.970     44.13    0.38
93 
+gold     0.0175    34.52    0.25
94 
+gold     0.500     42.74    0.20
95 
+gold     0.440     42.08    0.19
96 
+gold     0.0297    36.12    0.20
97 
+gold     0.0104    33.73    0.33
98 
+gold     0.778     43.81    0.35
99 
+gold     0.860     44.03    0.30
100 
+gold     0.538     42.66    0.18
101 
+gold     1.755     45.53    0.35
102 
+gold     0.886     42.91    0.81
103 
+gold     0.828     43.61    0.61
104 
+gold     0.630     42.62    0.24
105 
+gold     0.0170    34.47    0.23
106 
+gold     0.460     41.83    0.40
107 
+silver   0.638     43.30    0.36
108 
+gold     0.740     43.35    0.30
109 
+silver   0.644     42.78    0.26
110 
+gold     0.430     42.36    0.25
111 
+gold     0.0104    33.21    0.32
112 
+gold     0.0171    34.68    0.24
113 
+gold     0.0327    36.08    0.19
114 
+gold     0.053     36.97    0.18
115 
+gold     0.0170    34.18    0.23
116 
+gold     0.0234    35.36    0.20
117 
+gold     0.460     42.56    0.27
118 
+gold     0.500     42.75    0.19
119 
+gold     0.0257    35.41    0.20
120 
+gold     0.0157    34.58    0.24
121 
+gold     0.0316    35.85    0.19
122 
+gold     0.0104    33.56    0.31
123 
+silver   0.0121    34.05    0.32
124 
+gold     0.0141    34.43    0.26
125 
+gold     0.0136    33.73    0.26
126 
+gold     0.0380    36.67    0.18
127 
+gold     0.278     41.00    0.41
128 
+gold     1.056     44.25    0.23
129 
+gold     0.949     43.99    0.25
130 
+gold     0.815     43.76    0.33
131 
+gold     0.455     42.29    0.28
132 
+gold     1.19	   44.19    0.34
133 
+silver   0.369     41.62    0.31
134 
+gold     0.477     42.38    0.21
135 
+gold     0.0260    35.62    0.20
136 
+gold     0.0193    34.59    0.23
137 
+gold     0.0266    35.36    0.21
138 
+gold     0.0360    36.39    0.18
139 
+gold     0.0233    35.14    0.21
140 
+silver   0.0164    34.47    0.23
141 
+gold     0.0164    34.41    0.24
142 
+gold     0.500     42.75    0.24
143 
+gold     0.490     42.41    0.25
144 
+gold     0.470     42.74    0.23
145 
+gold     0.540     41.96    0.41
146 
+silver   0.420     40.79    0.32
147 
+gold     0.470     42.77    0.21
148 
+gold     0.543     42.68    0.19
149 
+gold     0.0218    35.06    0.21
150 
+gold     0.0162    34.13    0.23
151 
+gold     0.873     43.75    0.38
152 
+gold     0.771     43.12    0.17
153 
+gold     0.811     43.97    0.35
154 
+gold     0.798     43.88    0.31
155 
+gold     0.832     43.55    0.29
156 
+gold     0.882     43.90    0.30
157 
+gold     0.340     40.71    0.27
158 
+gold     0.397     40.89    0.30
159 
+gold     0.570     42.88    0.31
160 
+gold     0.710     43.05    0.32
161 
+gold     0.528     42.77    0.25
162 
+gold     0.884     44.23    0.19
163 
+silver   0.698     43.33    0.32
164 
+gold     0.815     44.09    0.28
165 
+gold     0.977     43.91    0.26
166 
+silver   0.935     43.99    0.38
167 
+silver   0.719     43.22    0.26
168 
+silver   0.422     42.02    0.17
169 
+silver   0.514     42.39    0.27
170 
+gold     0.475     42.14    0.19
171 
+gold     0.95	   44.06    0.26
172 
+gold     1.30	   45.27    0.19
173 
+silver   1.40	   45.09    0.45
174 
+gold     0.526     43.01    0.27
175 
+gold     1.305     44.70    0.22
176 
+silver   0.216     40.33    0.18
177 
+gold     0.735     43.09    0.19
178 
+gold     1.140     44.84    0.30
179 
+gold     1.265     45.20    0.20
180 
+gold     1.551     45.30    0.22
181 
+gold     0.67	   43.19    0.28
182 
+gold     0.64	   43.07    0.21
183 
+gold     1.340     45.05    0.25
184 
+gold     0.954     44.28    0.31
185 
+gold     0.839     43.86    0.22
186 
+gold     0.899     43.64    0.25
187 
+gold     0.94	   43.87    0.20
6/solution_sheet_6.R
Zeige Datei @ 052b6b9
... ... 
@@ -0,0 +1,236 @@
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
233 
+pOmarg=function(Oi) {I=integrate(posterior,13.0,18.0,Omegam=Oi)
234 
+    return(I$value)
235 
+    }
236 
+
0 237