finaly (8b82f20) - SS2012-StatisticalMethods.git

7/solutions7.r

@@ -1,8 +1,9 @@
 #Problem 1
-
+library(Hmisc)
 #definition of a straight line function
 linrel=function(x,m,b) m*x+b
-
+#sigvar = 0.4
+sigvar = 1
 #definition x-Vctor
 x_vec=seq(-1,3,0.1)
 
@@ -15,7 +16,7 @@ plot(x_vec, y_vec, type="l", col="black", xlab='x',ylab='f(x)=y')
 equ_points=10
 mock_vec = rep(NA,10)
 x_vec_2 =seq(0,2, 2/(equ_points-1))
-mock_vec = rnorm({1:10},mean = x_vec_2, sd =  0.4)
+mock_vec = rnorm({1:10},mean = x_vec_2, sd =  sigvar)
 points(x_vec_2,mock_vec, col="red", type="p")
 
 ##Problem2
@@ -29,9 +30,9 @@ b_vec = seq(-2,2,length.out = reso)
 
 
 # Define function chi2_quad
-chi2_qua = function(a,m,b) sum((mock_vec - quarel(mock_vec,a,m,b))**2/(0.4)**2)
-chi2_lin = function(m,b) sum((mock_vec - linrel(mock_vec,m,b))**2/(0.4)**2)
-chi2_con = function(b) sum((mock_vec - conrel(b))**2/(0.4)**2)
+chi2_qua = function(a,m,b) sum((mock_vec - quarel(mock_vec,a,m,b))**2/(sigvar)**2)
+chi2_lin = function(m,b) sum((mock_vec - linrel(mock_vec,m,b))**2/(sigvar)**2)
+chi2_con = function(b) sum((mock_vec - conrel(b))**2/(sigvar)**2)
 
 # Calculate 3-dim grid for values a=-4..4, m=-4..8, b-2..2
 qua_arr = array(NA, dim=c(length(b_vec),length(m_vec),length(a_vec)))
@@ -90,42 +91,58 @@ err_lin_b = 0
 err_lin_m = 0
 err_con_b = 0
 
-#Marginalization of posterior_quad over a
+#Marginalization of posterior_quad over a and m
 int_qua_a=array(0, c(dim_m, dim_b))
 for (i in (1:dim_m)){
 	for (j in (1:dim_b)){
-		I=int(a_vec, norm(posterior_qua[j,i,]))
+		I=int(a_vec, (posterior_qua)[j,i,])
 		int_qua_a[i,j]=I
 	}
 }
 
-#Marginalization of posterior_quad over b
+int_qua_am=array(0, c(dim_b))
+for (i in (1:dim_b)){
+		I=int(m_vec, (int_qua_a)[,i])
+		int_qua_am[i]=I
+}
+#Marginalization of posterior_quad over b and m
 int_qua_b=array(0, c(dim_m, dim_a))
 for (i in (1:dim_m)){
 	for (j in (1:dim_a)){
-		I=int(b_vec, norm(posterior_qua[,i,j]))
+		I=int(b_vec, (posterior_qua)[,i,j])
 		int_qua_b[i,j]=I
 	}
 }
-#Marginalization of posterior_quad over m
-int_qua_m=array(0, c(dim_a, dim_b))
+
+int_qua_bm=array(0, c(dim_a))
 for (i in (1:dim_a)){
-	for (j in (1:dim_b)){
-		I=int(m_vec, norm(posterior_qua[j,,i]))
-		int_qua_m[i,j]=I
+		I=int(m_vec, (int_qua_b)[,i])
+		int_qua_bm[i]=I
 }
+
+#Marginalization of posterior_quad over b and a
+
+int_qua_ba=array(0, c(dim_m))
+for (i in (1:dim_m)){
+		I=int(b_vec, (int_qua_a)[,i])
+		int_qua_ba[i]=I
 }
 
+
+
+
+######
+
 #Marginalization of posterior_lin over b and m
 int_lin_b=array(0, dim=c(length(m_vec)))
 for (i in (1:length(m_vec))){
-	I=int(b_vec, norm(posterior_lin[i,]))
+	I=int(b_vec, (posterior_lin[i,]))
 		int_lin_b[i]=I
 }
 
 int_lin_m=array(0, dim=c(length(b_vec)))
 for (i in (1:length(b_vec))){
-	I=int(m_vec, norm(posterior_lin[,i]))
+	I=int(m_vec, (posterior_lin[,i]))
 		int_lin_m[i]=I
 }
 
@@ -133,18 +150,97 @@ for (i in (1:length(b_vec))){
 #Marginalization of posterior_con over b
 
 int_lin_m=array(0, dim=c(dim_b))
-I=int(b_vec, norm(posterior_con))
-int_con_b=I
+I=int(b_vec, (posterior_con))
+int_con_b=posterior_con
+
+#search for best sigma
+sdtest=seq(0.01,1,by=0.01)
+
+sigma_test=function(A,mu,sd,xvec,f){
+	chi=0
+	for(j in (1:length(xvec)))
+		{
+			chi=chi+(f[j]-(A*exp(-0.5*(((xvec[j]-mu)**2)/(sd**2))))/sigvar)**2
+        }
+        return(chi)
+    }
+
+# create grids for vaious sigma functions
+# later, find min of grid
+
+grid_con_b <-array(0, dim=c(length(sdtest)))
+grid_lin_b <-array(0, dim=c(length(sdtest)))
+grid_lin_m <-array(0, dim=c(length(sdtest)))
+grid_qua_b <-array(0, dim=c(length(sdtest)))
+grid_qua_m <-array(0, dim=c(length(sdtest)))
+grid_qua_a <-array(0, dim=c(length(sdtest)))
+
+for(i in (1:length(sdtest)))
+{
+	grid_con_b[i] = sigma_test(max(int_con_b),best_b_con,sdtest[i], b_vec,int_con_b)
+	grid_lin_b[i] = sigma_test(max(int_lin_b),best_b_lin,sdtest[i], b_vec,int_lin_b)
+	grid_lin_m[i] = sigma_test(max(int_lin_m),best_m_lin,sdtest[i], m_vec,int_lin_m)
+	grid_qua_b[i] = sigma_test(max(int_qua_am),best_b_qua,sdtest[i], b_vec,int_qua_am)
+	grid_qua_m[i] = sigma_test(max(int_qua_ba),best_m_qua,sdtest[i], m_vec,int_qua_ba)
+	grid_qua_a[i] = sigma_test(max(int_qua_bm),best_a_qua,sdtest[i], a_vec,int_qua_bm)
+}
+err_con_b=sdtest[which(grid_con_b==min(grid_con_b),arr.ind=TRUE)]
+err_lin_b=sdtest[which(grid_lin_b==min(grid_lin_b),arr.ind=TRUE)]
+err_lin_m=sdtest[which(grid_lin_m==min(grid_lin_m),arr.ind=TRUE)]
+err_qua_b=sdtest[which(grid_qua_b==min(grid_qua_b),arr.ind=TRUE)]
+err_qua_m=sdtest[which(grid_qua_m==min(grid_qua_m),arr.ind=TRUE)]
+err_qua_a=sdtest[which(grid_qua_a==min(grid_qua_a),arr.ind=TRUE)]
+err_lin_m = err_lin_m[1]
+
 
 
 
 # The integration via fit is computationally expensive compared to the simple numerical integration 
 
 
-strqua = paste("quadratic: ax^2 + mx +b with a=",toString(best_a_qua),"+-",toString(err_qua_a),",\n m=", toString(best_m_qua),"+-",toString(err_qua_m),", b=",toString(best_b_qua),"+-",toString(err_qua_b))
-strlin = paste("linear: mx +b with m=", toString(best_m_lin),"+-",toString(err_lin_m),",\n b=",toString(best_b_lin),"+-",toString(err_lin_b))
-strcon = paste("constant: b with b=",toString(best_b_con),"+-",toString(err_con_b))
+strqua = paste("quadratic: ax^2 + mx +b with a=",round(best_a_qua,2),"+-",round(err_qua_a,2),",\n m=", round(best_m_qua,2),"+-",round(err_qua_m,2),", b=",round(best_b_qua,2),"+-",round(err_qua_b,2))
+strlin = paste("linear: mx +b with m=", round(best_m_lin,2),"+-",round(err_lin_m,2),",\n b=",round(best_b_lin,2),"+-",round(err_lin_b,2))
+strcon = paste("constant: b with b=",round(best_b_con,2),"+-",round(err_con_b,2))
 print(strqua)
 text(0,2.8,strqua)
 text(0,2.1,strlin)
 text(0,1.6,strcon)
+
+#Plot errorbars quad
+for( x in (x_vec_2)){
+
+tmperrup = quarel(x,best_a_qua,best_m_qua,best_b_qua) + (err_qua_a*x**2+x*err_qua_m+err_qua_b)
+tmperrdo = quarel(x,best_a_qua,best_m_qua,best_b_qua) - (err_qua_a*x**2+x*err_qua_m+err_qua_b)
+errbar(x,quarel(x,best_a_qua,best_m_qua,best_b_qua),tmperrup,tmperrdo,add=TRUE)
+
+tmperrup = linrel(x,best_m_lin,best_b_lin) + (x*err_lin_m+err_lin_b)
+tmperrdo = linrel(x,best_m_lin,best_b_lin) - (x*err_lin_m+err_lin_b)
+errbar(x,linrel(x,best_m_lin,best_b_lin),tmperrup,tmperrdo,add=TRUE)
+
+tmperrup = best_b_lin + err_con_b
+tmperrdo = best_b_lin - err_con_b
+#print(tmperrup)
+#print(tmperrdo)
+#print(x)
+#print(best_b_con)
+errbar(x,best_b_con,tmperrup,tmperrdo,add=TRUE)
+}
+#Problem 4
+# calculate the evidence
+evidence_qua=err_qua_a*err_qua_m*err_qua_b*max(posterior_qua)*(1/sum(posterior_qua))
+evidence_lin=err_lin_m*err_lin_b*max(posterior_lin)*(1/sum(posterior_lin))
+evidence_con=err_con_b*max(posterior_con)* (1/sum(posterior_con))
+print('evidence')
+print('quadratic')
+print(evidence_qua)
+print('linear')
+print(evidence_lin)
+print('const.')
+print(evidence_con)
+
+#the constant would be the best model. This is in contrast to the naive expectation that straight line model would fit best.
+
+#Problem 5
+#To run the script with sigma 1 change the variable sigvar to 1
+#Theoretically the evidence should point to the quadratic model because it could better fit the new mock data with higher sigvar. We get the constant again as the best model. Forthermor we gat bigger error bars and so a better evidence.
+


...	...	@@ -1,8 +1,9 @@
1	1	#Problem 1
2		-
	2	+library(Hmisc)
3	3	#definition of a straight line function
4	4	linrel=function(x,m,b) m*x+b
5		-
	5	+#sigvar = 0.4
	6	+sigvar = 1
6	7	#definition x-Vctor
7	8	x_vec=seq(-1,3,0.1)
8	9
...	...	@@ -15,7 +16,7 @@ plot(x_vec, y_vec, type="l", col="black", xlab='x',ylab='f(x)=y')
15	16	equ_points=10
16	17	mock_vec = rep(NA,10)
17	18	x_vec_2 =seq(0,2, 2/(equ_points-1))
18		-mock_vec = rnorm({1:10},mean = x_vec_2, sd = 0.4)
	19	+mock_vec = rnorm({1:10},mean = x_vec_2, sd = sigvar)
19	20	points(x_vec_2,mock_vec, col="red", type="p")
20	21
21	22	##Problem2
...	...	@@ -29,9 +30,9 @@ b_vec = seq(-2,2,length.out = reso)
29	30
30	31
31	32	# Define function chi2_quad
32		-chi2_qua = function(a,m,b) sum((mock_vec - quarel(mock_vec,a,m,b))2/(0.4)2)
33		-chi2_lin = function(m,b) sum((mock_vec - linrel(mock_vec,m,b))2/(0.4)2)
34		-chi2_con = function(b) sum((mock_vec - conrel(b))2/(0.4)2)
	33	+chi2_qua = function(a,m,b) sum((mock_vec - quarel(mock_vec,a,m,b))2/(sigvar)2)
	34	+chi2_lin = function(m,b) sum((mock_vec - linrel(mock_vec,m,b))2/(sigvar)2)
	35	+chi2_con = function(b) sum((mock_vec - conrel(b))2/(sigvar)2)
35	36
36	37	# Calculate 3-dim grid for values a=-4..4, m=-4..8, b-2..2
37	38	qua_arr = array(NA, dim=c(length(b_vec),length(m_vec),length(a_vec)))
...	...	@@ -90,42 +91,58 @@ err_lin_b = 0
90	91	err_lin_m = 0
91	92	err_con_b = 0
92	93
93		-#Marginalization of posterior_quad over a
	94	+#Marginalization of posterior_quad over a and m
94	95	int_qua_a=array(0, c(dim_m, dim_b))
95	96	for (i in (1:dim_m)){
96	97	for (j in (1:dim_b)){
97		- I=int(a_vec, norm(posterior_qua[j,i,]))
	98	+ I=int(a_vec, (posterior_qua)[j,i,])
98	99	int_qua_a[i,j]=I
99	100	}
100	101	}
101	102
102		-#Marginalization of posterior_quad over b
	103	+int_qua_am=array(0, c(dim_b))
	104	+for (i in (1:dim_b)){
	105	+ I=int(m_vec, (int_qua_a)[,i])
	106	+ int_qua_am[i]=I
	107	+}
	108	+#Marginalization of posterior_quad over b and m
103	109	int_qua_b=array(0, c(dim_m, dim_a))
104	110	for (i in (1:dim_m)){
105	111	for (j in (1:dim_a)){
106		- I=int(b_vec, norm(posterior_qua[,i,j]))
	112	+ I=int(b_vec, (posterior_qua)[,i,j])
107	113	int_qua_b[i,j]=I
108	114	}
109	115	}
110		-#Marginalization of posterior_quad over m
111		-int_qua_m=array(0, c(dim_a, dim_b))
	116	+
	117	+int_qua_bm=array(0, c(dim_a))
112	118	for (i in (1:dim_a)){
113		- for (j in (1:dim_b)){
114		- I=int(m_vec, norm(posterior_qua[j,,i]))
115		- int_qua_m[i,j]=I
	119	+ I=int(m_vec, (int_qua_b)[,i])
	120	+ int_qua_bm[i]=I
116	121	}
	122	+
	123	+#Marginalization of posterior_quad over b and a
	124	+
	125	+int_qua_ba=array(0, c(dim_m))
	126	+for (i in (1:dim_m)){
	127	+ I=int(b_vec, (int_qua_a)[,i])
	128	+ int_qua_ba[i]=I
117	129	}
118	130
	131	+
	132	+
	133	+
	134	+######
	135	+
119	136	#Marginalization of posterior_lin over b and m
120	137	int_lin_b=array(0, dim=c(length(m_vec)))
121	138	for (i in (1:length(m_vec))){
122		- I=int(b_vec, norm(posterior_lin[i,]))
	139	+ I=int(b_vec, (posterior_lin[i,]))
123	140	int_lin_b[i]=I
124	141	}
125	142
126	143	int_lin_m=array(0, dim=c(length(b_vec)))
127	144	for (i in (1:length(b_vec))){
128		- I=int(m_vec, norm(posterior_lin[,i]))
	145	+ I=int(m_vec, (posterior_lin[,i]))
129	146	int_lin_m[i]=I
130	147	}
131	148
...	...	@@ -133,18 +150,97 @@ for (i in (1:length(b_vec))){
133	150	#Marginalization of posterior_con over b
134	151
135	152	int_lin_m=array(0, dim=c(dim_b))
136		-I=int(b_vec, norm(posterior_con))
137		-int_con_b=I
	153	+I=int(b_vec, (posterior_con))
	154	+int_con_b=posterior_con
	155	+
	156	+#search for best sigma
	157	+sdtest=seq(0.01,1,by=0.01)
	158	+
	159	+sigma_test=function(A,mu,sd,xvec,f){
	160	+ chi=0
	161	+ for(j in (1:length(xvec)))
	162	+ {
	163	+ chi=chi+(f[j]-(Aexp(-0.5(((xvec[j]-mu)2)/(sd2))))/sigvar)**2
	164	+ }
	165	+ return(chi)
	166	+ }
	167	+
	168	+# create grids for vaious sigma functions
	169	+# later, find min of grid
	170	+
	171	+grid_con_b <-array(0, dim=c(length(sdtest)))
	172	+grid_lin_b <-array(0, dim=c(length(sdtest)))
	173	+grid_lin_m <-array(0, dim=c(length(sdtest)))
	174	+grid_qua_b <-array(0, dim=c(length(sdtest)))
	175	+grid_qua_m <-array(0, dim=c(length(sdtest)))
	176	+grid_qua_a <-array(0, dim=c(length(sdtest)))
	177	+
	178	+for(i in (1:length(sdtest)))
	179	+{
	180	+ grid_con_b[i] = sigma_test(max(int_con_b),best_b_con,sdtest[i], b_vec,int_con_b)
	181	+ grid_lin_b[i] = sigma_test(max(int_lin_b),best_b_lin,sdtest[i], b_vec,int_lin_b)
	182	+ grid_lin_m[i] = sigma_test(max(int_lin_m),best_m_lin,sdtest[i], m_vec,int_lin_m)
	183	+ grid_qua_b[i] = sigma_test(max(int_qua_am),best_b_qua,sdtest[i], b_vec,int_qua_am)
	184	+ grid_qua_m[i] = sigma_test(max(int_qua_ba),best_m_qua,sdtest[i], m_vec,int_qua_ba)
	185	+ grid_qua_a[i] = sigma_test(max(int_qua_bm),best_a_qua,sdtest[i], a_vec,int_qua_bm)
	186	+}
	187	+err_con_b=sdtest[which(grid_con_b==min(grid_con_b),arr.ind=TRUE)]
	188	+err_lin_b=sdtest[which(grid_lin_b==min(grid_lin_b),arr.ind=TRUE)]
	189	+err_lin_m=sdtest[which(grid_lin_m==min(grid_lin_m),arr.ind=TRUE)]
	190	+err_qua_b=sdtest[which(grid_qua_b==min(grid_qua_b),arr.ind=TRUE)]
	191	+err_qua_m=sdtest[which(grid_qua_m==min(grid_qua_m),arr.ind=TRUE)]
	192	+err_qua_a=sdtest[which(grid_qua_a==min(grid_qua_a),arr.ind=TRUE)]
	193	+err_lin_m = err_lin_m[1]
	194	+
138	195
139	196
140	197
141	198	# The integration via fit is computationally expensive compared to the simple numerical integration
142	199
143	200
144		-strqua = paste("quadratic: ax^2 + mx +b with a=",toString(best_a_qua),"+-",toString(err_qua_a),",\n m=", toString(best_m_qua),"+-",toString(err_qua_m),", b=",toString(best_b_qua),"+-",toString(err_qua_b))
145		-strlin = paste("linear: mx +b with m=", toString(best_m_lin),"+-",toString(err_lin_m),",\n b=",toString(best_b_lin),"+-",toString(err_lin_b))
146		-strcon = paste("constant: b with b=",toString(best_b_con),"+-",toString(err_con_b))
	201	+strqua = paste("quadratic: ax^2 + mx +b with a=",round(best_a_qua,2),"+-",round(err_qua_a,2),",\n m=", round(best_m_qua,2),"+-",round(err_qua_m,2),", b=",round(best_b_qua,2),"+-",round(err_qua_b,2))
	202	+strlin = paste("linear: mx +b with m=", round(best_m_lin,2),"+-",round(err_lin_m,2),",\n b=",round(best_b_lin,2),"+-",round(err_lin_b,2))
	203	+strcon = paste("constant: b with b=",round(best_b_con,2),"+-",round(err_con_b,2))
147	204	print(strqua)
148	205	text(0,2.8,strqua)
149	206	text(0,2.1,strlin)
150	207	text(0,1.6,strcon)
	208	+
	209	+#Plot errorbars quad
	210	+for( x in (x_vec_2)){
	211	+
	212	+tmperrup = quarel(x,best_a_qua,best_m_qua,best_b_qua) + (err_qua_ax2+xerr_qua_m+err_qua_b)
	213	+tmperrdo = quarel(x,best_a_qua,best_m_qua,best_b_qua) - (err_qua_ax2+xerr_qua_m+err_qua_b)
	214	+errbar(x,quarel(x,best_a_qua,best_m_qua,best_b_qua),tmperrup,tmperrdo,add=TRUE)
	215	+
	216	+tmperrup = linrel(x,best_m_lin,best_b_lin) + (x*err_lin_m+err_lin_b)
	217	+tmperrdo = linrel(x,best_m_lin,best_b_lin) - (x*err_lin_m+err_lin_b)
	218	+errbar(x,linrel(x,best_m_lin,best_b_lin),tmperrup,tmperrdo,add=TRUE)
	219	+
	220	+tmperrup = best_b_lin + err_con_b
	221	+tmperrdo = best_b_lin - err_con_b
	222	+#print(tmperrup)
	223	+#print(tmperrdo)
	224	+#print(x)
	225	+#print(best_b_con)
	226	+errbar(x,best_b_con,tmperrup,tmperrdo,add=TRUE)
	227	+}
	228	+#Problem 4
	229	+# calculate the evidence
	230	+evidence_qua=err_qua_aerr_qua_merr_qua_bmax(posterior_qua)(1/sum(posterior_qua))
	231	+evidence_lin=err_lin_merr_lin_bmax(posterior_lin)*(1/sum(posterior_lin))
	232	+evidence_con=err_con_bmax(posterior_con) (1/sum(posterior_con))
	233	+print('evidence')
	234	+print('quadratic')
	235	+print(evidence_qua)
	236	+print('linear')
	237	+print(evidence_lin)
	238	+print('const.')
	239	+print(evidence_con)
	240	+
	241	+#the constant would be the best model. This is in contrast to the naive expectation that straight line model would fit best.
	242	+
	243	+#Problem 5
	244	+#To run the script with sigma 1 change the variable sigvar to 1
	245	+#Theoretically the evidence should point to the quadratic model because it could better fit the new mock data with higher sigvar. We get the constant again as the best model. Forthermor we gat bigger error bars and so a better evidence.
	246	+
151	247