SS2012-StatisticalMethods.git

1) #Problem 1
2) 
3) #definition of a straight line function

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

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

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

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

16) mock_vec = rep(NA,10)
17) x_vec_2 =seq(0,2, 2/(equ_points-1))
18) mock_vec = rnorm({1:10},mean = x_vec_2, sd =  0.4)
19) points(x_vec_2,mock_vec, col="red", type="p")
20) 
21) ##Problem2
22) conrel =function(b) b
23) quarel =function(x,a,m,b) a*x**2 + m*x+b
24) 

25) reso = 50  #the resolution

26) a_vec = seq(-4,4,length.out = reso)
27) m_vec = seq(-4,8,length.out = reso)
28) b_vec = seq(-2,2,length.out = reso)
29) 
30) 
31) # 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)
35) 
36) # Calculate 3-dim grid for values a=-4..4, m=-4..8, b-2..2
37) qua_arr = array(NA, dim=c(length(b_vec),length(m_vec),length(a_vec)))
38) lin_arr = array(NA, dim=c(length(b_vec),length(m_vec)))
39) con_arr = array(NA, dim=c(length(b_vec)))
40) for (i in (1:length(b_vec))){
41)     con_arr[i] = chi2_con(b_vec[i])
42) 	for (j in (1:length(m_vec))){
43)         lin_arr[i,j] = chi2_lin(m_vec[j],b_vec[i])
44) 		for (k in (1:length(a_vec))){
45) 			qua_arr[i,j,k] = chi2_qua(a_vec[k],m_vec[j],b_vec[i])
46) 		}
47) 	}
48) }

49) #find min

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

54) 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')
55) 

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

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

60) minimum_con=min(con_arr)

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

62) points(x_vec, rep(conrel(best_b_lin),41), type="l", col="yellow", xlab='x',ylab='f(x)=y')
63) 
64) #Probelm3
65) norm = function(a) 1/sum(a) * a
66) 

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

70) 

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

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

93) #Marginalization of posterior_quad over a
94) int_qua_a=array(0, c(dim_m, dim_b))
95) for (i in (1:dim_m)){
96) 	for (j in (1:dim_b)){
97) 		I=int(a_vec, norm(posterior_qua[j,i,]))
98) 		int_qua_a[i,j]=I
99) 	}
100) }
101) 
102) #Marginalization of posterior_quad over b
103) int_qua_b=array(0, c(dim_m, dim_a))
104) for (i in (1:dim_m)){
105) 	for (j in (1:dim_a)){
106) 		I=int(b_vec, norm(posterior_qua[,i,j]))
107) 		int_qua_b[i,j]=I
108) 	}
109) }
110) #Marginalization of posterior_quad over m
111) int_qua_m=array(0, c(dim_a, dim_b))
112) 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
116) 	}
117) }
118) 
119) #Marginalization of posterior_lin over b and m
120) int_lin_b=array(0, dim=c(length(m_vec)))
121) for (i in (1:length(m_vec))){
122) 	I=int(b_vec, norm(posterior_lin[i,]))
123) 		int_lin_b[i]=I
124) }
125) 
126) int_lin_m=array(0, dim=c(length(b_vec)))
127) for (i in (1:length(b_vec))){
128) 	I=int(m_vec, norm(posterior_lin[,i]))
129) 		int_lin_m[i]=I
130) }
131) 
132) 
133) #Marginalization of posterior_con over b
134) 
135) int_lin_m=array(0, dim=c(dim_b))
136) I=int(b_vec, norm(posterior_con))
137) int_con_b=I