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(lin_arr)
61) best_b_con = b_vec[which(con_arr == minimum_con,arr.ind=TRUE)[1]]