8b82f20d5ded236bd5c65a48a39214bf190b0e22
 add sheet 7 Michi authored 12 years ago ``````1) #Problem 1 `````` finaly Michi authored 12 years ago ``````2) library(Hmisc) `````` add sheet 7 Michi authored 12 years ago ``````3) #definition of a straight line function `````` add 2.2 Michi authored 12 years ago ``````4) linrel=function(x,m,b) m*x+b `````` finaly Michi authored 12 years ago ``````5) #sigvar = 0.4 6) sigvar = 1 `````` add sheet 7 Michi authored 12 years ago ``````7) #definition x-Vctor 8) x_vec=seq(-1,3,0.1) 9) 10) #definition y-vector `````` add 2.2 Michi authored 12 years ago ``````11) y_vec=linrel(x_vec,1,0) `````` add sheet 7 Michi authored 12 years ago ``````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 `````` new 2 i Michi authored 12 years ago ``````17) mock_vec = rep(NA,10) 18) x_vec_2 =seq(0,2, 2/(equ_points-1)) `````` finaly Michi authored 12 years ago ``````19) mock_vec = rnorm({1:10},mean = x_vec_2, sd = sigvar) `````` new 2 i Michi authored 12 years ago ``````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) `````` to 2.3 Michi authored 12 years ago ``````26) reso = 50 #the resolution `````` add 2.2 Michi authored 12 years ago ``````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 `````` finaly Michi authored 12 years ago ``````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) `````` add 2.2 Michi authored 12 years ago ``````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) } `````` to 2.3 Michi authored 12 years ago ``````50) #find min `````` add 2.2 Michi authored 12 years ago ``````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]] `````` to 2.3 Michi authored 12 years ago ``````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) `````` add 2.2 Michi authored 12 years ago ``````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]] `````` to 2.3 Michi authored 12 years ago ``````60) points(x_vec, linrel(x_vec,best_m_lin,best_b_lin), type="l", col="blue", xlab='x',ylab='f(x)=y') `````` add marginalized Michi authored 12 years ago ``````61) minimum_con=min(con_arr) `````` add 2.2 Michi authored 12 years ago ``````62) best_b_con = b_vec[which(con_arr == minimum_con,arr.ind=TRUE)[1]] `````` to 2.3 Michi authored 12 years ago ``````63) points(x_vec, rep(conrel(best_b_lin),41), type="l", col="yellow", xlab='x',ylab='f(x)=y') 64) 65) #Probelm3 66) norm = function(a) 1/sum(a) * a 67) `````` add marginalized Michi authored 12 years ago ``````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)])) `````` to 2.3 Michi authored 12 years ago ``````71) `````` add marginalized Michi authored 12 years ago ``````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) `````` new 2.3 Michi authored 12 years ago ``````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) `````` finaly Michi authored 12 years ago ``````94) #Marginalization of posterior_quad over a and m `````` add marginalized Michi authored 12 years ago ``````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)){ `````` finaly Michi authored 12 years ago ``````98) I=int(a_vec, (posterior_qua)[j,i,]) `````` add marginalized Michi authored 12 years ago ``````99) int_qua_a[i,j]=I 100) } 101) } 102) `````` finaly Michi authored 12 years ago ``````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 `````` add marginalized Michi authored 12 years ago ``````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)){ `````` finaly Michi authored 12 years ago ``````112) I=int(b_vec, (posterior_qua)[,i,j]) `````` add marginalized Michi authored 12 years ago ``````113) int_qua_b[i,j]=I 114) } 115) } `````` finaly Michi authored 12 years ago ``````116) 117) int_qua_bm=array(0, c(dim_a)) `````` add marginalized Michi authored 12 years ago ``````118) for (i in (1:dim_a)){ `````` finaly Michi authored 12 years ago ``````119) I=int(m_vec, (int_qua_b)[,i]) 120) int_qua_bm[i]=I `````` add marginalized Michi authored 12 years ago ``````121) } 122) `````` finaly Michi authored 12 years ago ``````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) `````` add marginalized Michi authored 12 years ago ``````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))){ `````` finaly Michi authored 12 years ago ``````139) I=int(b_vec, (posterior_lin[i,])) `````` add marginalized Michi authored 12 years ago ``````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))){ `````` finaly Michi authored 12 years ago ``````145) I=int(m_vec, (posterior_lin[,i])) `````` add marginalized Michi authored 12 years ago ``````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)) `````` finaly Michi authored 12 years ago ``````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) `````` new 2.3 Michi authored 12 years ago ``````195) 196) 197) 198) # The integration via fit is computationally expensive compared to the simple numerical integration 199) 200) `````` finaly Michi authored 12 years ago ``````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)) `````` new 2.3 Michi authored 12 years ago ``````204) print(strqua) 205) text(0,2.8,strqua) 206) text(0,2.1,strlin) 207) text(0,1.6,strcon) ``````