SS2012-StatisticalMethods.git

The fist commit
Michi

Michi commited on 2012-05-14 14:57:17
Zeige 4 geänderte Dateien mit 319 Einfügungen und 0 Löschungen.

4/coursework4.pdf
Zeige Datei @ a7fc161
4/dL.R
Zeige Datei @ a7fc161
... ... 
@@ -0,0 +1,25 @@
1 
+#Constants
2 
+H0=72.0 #km/s/Mpc
3 
+c = 3*10^5 #km/s
4 
+
5 
+
6 
+#Hubble's law
7 
+Hubble=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))
8 
+
9 
+
10 
+#H^(-1)
11 
+Hinv=function(zint,Omegam) 1.0/Hubble(zint,Omegam)
12 
+
13 
+
14 
+#Luminosity distance
15 
+dL=function(z,Omegami){
16 
+	dLsol=z
17 
+	for (i in (1:length(z))){
18 
+		zarg=z[i]
19 
+		I=integrate(Hinv,0.0,zarg,Omegam=Omegami)
20 
+		dLsol[i] = c*(1+zarg)*I$value
21 
+		}
22 
+	return(dLsol)
23 
+	}
24 
+
25 
+
4/sn_data_riess.dat
Zeige Datei @ a7fc161
... ... 
@@ -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
4/soluton.r
Zeige Datei @ a7fc161
... ... 
@@ -0,0 +1,107 @@
1 
+#Exercise 1
2 
+#a)
3 
+#setwd("C:/Users/Studium/Desktop/Statistische Methoden/Sheet 4")
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 
+data
7 
+
8 
+#b)
9 
+#Constants
10 
+H0=72.0 #km/s/Mpc
11 
+c = 3*10^5 #km/s
12 
+
13 
+
14 
+#Hubble's law
15 
+Hubble=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))
16 
+
17 
+
18 
+#H^(-1)
19 
+Hinv=function(zint,Omegam) 1.0/Hubble(zint,Omegam)
20 
+
21 
+
22 
+#Luminosity distance
23 
+dL=function(z,Omegami){
24 
+	dLsol=z
25 
+	for (i in (1:length(z))){
26 
+		zarg=z[i]
27 
+		I=integrate(Hinv,0.0,zarg,Omegam=Omegami)
28 
+		dLsol[i] = c*(1+zarg)*I$value
29 
+		}
30 
+	return(dLsol)
31 
+	}
32 
+
33 
+m=function(z, Omegam, M) M + 5*log10(H0 * dL(z, Omegam))
34 
+
35 
+mag_th=NULL
36 
+mag_th[186]=1
37 
+for (i in (1:186)){
38 
+	mag_th[i] = m(data[i,1], 1, data[i,2]) #Omegam = 1 for a flat universe
39 
+	}
40 
+
41 
+#c)
42 
+chisquare=function(Omegam, M) {
43 
+	#summand = (data[i,2] - m(data[i,1], Omegam, M))**2/data[i,3]**2
44 
+	chisqu=0
45 
+	for (i in (1:186)){
46 
+		chisqu=chisqu+(data[i,2] - m(data[i,1], Omegam, M))**2/data[i,3]**2
47 
+	}
48 
+	return(chisqu)
49 
+}
50 
+
51 
+
52 
+#d)
53 
+#Erzeugen der zwei Vektoren f�r schrittweises Wachsen von Omegam und M
54 
+Omegam_vec=seq(0.0,1.0, by=0.05)
55 
+M_vec=seq(15.5, 16.6, 0.01)
56 
+Omegam_vec
57 
+M_vec
58 
+
59 
+#Startwerte
60 
+chi_min=999999999
61 
+Omegam_min=0
62 
+M_min=0
63 
+
64 
+#zweifach verschachtelte Schleife f�r Omegam und M
65 
+chisquare_min = chisquare(0,0) #set initial value
66 
+for (i in (1: length(M_vec))){
67 
+	for (j in Omegam_vec){
68 
+		if (chisquare(j,i) < chisquare_min){ #compare current value with known smallest value
69 
+			chisquare_min = chisquare(j,i)# set new smallest value
70 
+			Omegam_min = j
71 
+			M_min = i
72 
+		}
73 
+	}
74 
+	#chisquare_min
75 
+	#Omegam_min
76 
+	#M_min
77 
+}
78 
+
79 
+#e)
80 
+posterior = function(Omegam, M){
81 
+	p = exp(-0.5*(chisquare(Omegam, M)- chisquare_min))
82 
+	return(p)
83 
+}
84 
+
85 
+#f)
86 
+
87 
+
88 
+jmax=length(Omegam_vec)
89 
+kmax=length(M_vec)
90 
+j=1
91 
+k=1
92 
+# generate a posterior matrix for the contour plot
93 
+postplot = matrix(nrow=jmax,ncol=kmax)
94 
+while (j<=jmax){ k=1
95 
+	while (k<=kmax) {
96 
+	postplot[j,k]=posterior(Omegam_vec[j],M_vec[k])
97 
+	k=k+1}
98 
+	j=j+1}
99 
+quartz()
100 
+contour(Omegam_vec,M_vec,postplot,drawlabels=FALSE,xlab='Omega',ylab='M')
101 
+
102 
+
103 
+Mmarg=function(Ai) {I=integrate(posterior,13,18,A=Ai)
104 
+	return(I$value)
105 
+	}
106 
+
107 
+
0 108