SS2012-StatisticalMethods.git

1) #Problem 1

2) library(Hmisc)

3) #definition of a straight line function

4) linrel=function(x,m,b) m*x+b

5) #sigvar = 0.4
6) sigvar = 1

7) #definition x-Vctor
8) x_vec=seq(-1,3,0.1)
9) 
10) #definition y-vector

11) y_vec=linrel(x_vec,1,0)

12) #plotting x-vector vs. y-vector
13) plot(x_vec, y_vec, type="l", col="black", xlab='x',ylab='f(x)=y')
14) 
15) #generating mock data point for 0 <= x <= 2 for 10 equaly spaced points
16) equ_points=10

17) mock_vec = rep(NA,10)
18) x_vec_2 =seq(0,2, 2/(equ_points-1))

19) mock_vec = rnorm({1:10},mean = x_vec_2, sd =  sigvar)

20) points(x_vec_2,mock_vec, col="red", type="p")
21) 
22) ##Problem2
23) conrel =function(b) b
24) quarel =function(x,a,m,b) a*x**2 + m*x+b
25) 

26) reso = 50  #the resolution

27) a_vec = seq(-4,4,length.out = reso)
28) m_vec = seq(-4,8,length.out = reso)
29) b_vec = seq(-2,2,length.out = reso)
30) 
31) 
32) # Define function chi2_quad

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)

36) 
37) # Calculate 3-dim grid for values a=-4..4, m=-4..8, b-2..2
38) qua_arr = array(NA, dim=c(length(b_vec),length(m_vec),length(a_vec)))
39) lin_arr = array(NA, dim=c(length(b_vec),length(m_vec)))
40) con_arr = array(NA, dim=c(length(b_vec)))
41) for (i in (1:length(b_vec))){
42)     con_arr[i] = chi2_con(b_vec[i])
43) 	for (j in (1:length(m_vec))){
44)         lin_arr[i,j] = chi2_lin(m_vec[j],b_vec[i])
45) 		for (k in (1:length(a_vec))){
46) 			qua_arr[i,j,k] = chi2_qua(a_vec[k],m_vec[j],b_vec[i])
47) 		}
48) 	}
49) }

50) #find min

51) minimum_qua=min(qua_arr)
52) best_m_qua = m_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[2]]
53) best_b_qua = b_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[1]]
54) best_a_qua = a_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[3]]

55) points(x_vec, quarel(x_vec,best_a_qua,best_m_qua,best_b_qua), type="l", col="red", xlab='x',ylab='f(x)=y')
56) 

57) minimum_lin=min(lin_arr)
58) best_m_lin = m_vec[which(lin_arr == minimum_lin,arr.ind=TRUE)[2]]
59) best_b_lin = b_vec[which(lin_arr == minimum_lin,arr.ind=TRUE)[1]]

60) points(x_vec, linrel(x_vec,best_m_lin,best_b_lin), type="l", col="blue", xlab='x',ylab='f(x)=y')

61) minimum_con=min(con_arr)

62) best_b_con = b_vec[which(con_arr == minimum_con,arr.ind=TRUE)[1]]

63) points(x_vec, rep(conrel(best_b_con),41), type="l", col="yellow", xlab='x',ylab='f(x)=y')

64) 
65) #Probelm3
66) norm = function(a) 1/sum(a) * a
67) 

68) posterior_qua =  exp(-0.5 *(norm(qua_arr) - norm(qua_arr)[which(qua_arr == minimum_qua,arr.ind=TRUE)]))
69) posterior_lin =  exp(-0.5 *(norm(lin_arr) - norm(lin_arr)[which(lin_arr == minimum_lin,arr.ind=TRUE)]))
70) posterior_con =  exp(-0.5 *(norm(con_arr) - norm(con_arr)[which(con_arr == minimum_con,arr.ind=TRUE)]))

71) 

72) 
73) #integration numerically function
74) int = function(x_vec, y_vec){
75) 	A=0
76) 	step=(max(x_vec)-min(x_vec))/(length(x_vec))
77) 	for (i in y_vec){
78) 		A=A+step*i
79) 		}
80) 	return(A)
81) }
82) 
83) dim_m = length(m_vec)
84) dim_a = length(a_vec)
85) dim_b = length(b_vec)
86) 

87) err_qua_m = 0
88) err_qua_a = 0
89) err_qua_b = 0 
90) err_lin_b = 0
91) err_lin_m = 0
92) err_con_b = 0
93) 

94) #Marginalization of posterior_quad over a and m

95) int_qua_a=array(0, c(dim_m, dim_b))
96) for (i in (1:dim_m)){
97) 	for (j in (1:dim_b)){

98) 		I=int(a_vec, (posterior_qua)[j,i,])

99) 		int_qua_a[i,j]=I
100) 	}
101) }
102) 

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

109) int_qua_b=array(0, c(dim_m, dim_a))
110) for (i in (1:dim_m)){
111) 	for (j in (1:dim_a)){

112) 		I=int(b_vec, (posterior_qua)[,i,j])

113) 		int_qua_b[i,j]=I
114) 	}
115) }

116) 
117) int_qua_bm=array(0, c(dim_a))

118) for (i in (1:dim_a)){

119) 		I=int(m_vec, (int_qua_b)[,i])
120) 		int_qua_bm[i]=I

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
129) }
130) 
131) 
132) 
133) 
134) ######
135) 

136) #Marginalization of posterior_lin over b and m
137) int_lin_b=array(0, dim=c(length(m_vec)))
138) for (i in (1:length(m_vec))){

139) 	I=int(b_vec, (posterior_lin[i,]))

140) 		int_lin_b[i]=I
141) }
142) 
143) int_lin_m=array(0, dim=c(length(b_vec)))
144) for (i in (1:length(b_vec))){

145) 	I=int(m_vec, (posterior_lin[,i]))

146) 		int_lin_m[i]=I
147) }
148) 
149) 
150) #Marginalization of posterior_con over b
151) 
152) int_lin_m=array(0, dim=c(dim_b))

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]-(A*exp(-0.5*(((xvec[j]-mu)**2)/(sd**2))))/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) 

195) 
196) 
197) 

198) # The integration via fit (spinefun) is computationally expensive compared to the simple numerical integration by summing up. 

199) 
200) 

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))

204) print(strqua)
205) text(0,2.8,strqua)
206) text(0,2.1,strlin)
207) text(0,1.6,strcon)

208) 

209) #Plot errorbars quad. Here i'm not sure how to plot the error bars so I plot them like done below. An explanation about how to plot error bar in the tutorial would be nice.