7/solutions7.r (90df4f2) - SS2012-StatisticalMethods.git

solutions7.r

#Problem 1

#definition of a straight line function
linrel=function(x,m,b) m*x+b

#definition x-Vctor
x_vec=seq(-1,3,0.1)

#definition y-vector
y_vec=linrel(x_vec,1,0)
#plotting x-vector vs. y-vector
plot(x_vec, y_vec, type="l", col="black", xlab='x',ylab='f(x)=y')

#generating mock data point for 0 <= x <= 2 for 10 equaly spaced points
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)
points(x_vec_2,mock_vec, col="red", type="p")

##Problem2
conrel =function(b) b
quarel =function(x,a,m,b) a*x**2 + m*x+b

reso = 50  #the resolution
a_vec = seq(-4,4,length.out = reso)
m_vec = seq(-4,8,length.out = reso)
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)

# 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)))
lin_arr = array(NA, dim=c(length(b_vec),length(m_vec)))
con_arr = array(NA, dim=c(length(b_vec)))
for (i in (1:length(b_vec))){
    con_arr[i] = chi2_con(b_vec[i])
	for (j in (1:length(m_vec))){
        lin_arr[i,j] = chi2_lin(m_vec[j],b_vec[i])
		for (k in (1:length(a_vec))){
			qua_arr[i,j,k] = chi2_qua(a_vec[k],m_vec[j],b_vec[i])
		}
	}
}
#find min
minimum_qua=min(qua_arr)
best_m_qua = m_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[2]]
best_b_qua = b_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[1]]
best_a_qua = a_vec[which(qua_arr == minimum_qua,arr.ind=TRUE)[3]]
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')

minimum_lin=min(lin_arr)
best_m_lin = m_vec[which(lin_arr == minimum_lin,arr.ind=TRUE)[2]]
best_b_lin = b_vec[which(lin_arr == minimum_lin,arr.ind=TRUE)[1]]
points(x_vec, linrel(x_vec,best_m_lin,best_b_lin), type="l", col="blue", xlab='x',ylab='f(x)=y')
minimum_con=min(con_arr)
best_b_con = b_vec[which(con_arr == minimum_con,arr.ind=TRUE)[1]]
points(x_vec, rep(conrel(best_b_lin),41), type="l", col="yellow", xlab='x',ylab='f(x)=y')

#Probelm3
norm = function(a) 1/sum(a) * a

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

#integration numerically function
int = function(x_vec, y_vec){
	A=0
	step=(max(x_vec)-min(x_vec))/(length(x_vec))
	for (i in y_vec){
		A=A+step*i
		}
	return(A)
}

dim_m = length(m_vec)
dim_a = length(a_vec)
dim_b = length(b_vec)

err_qua_m = 0
err_qua_a = 0
err_qua_b = 0 
err_lin_b = 0
err_lin_m = 0
err_con_b = 0

#Marginalization of posterior_quad over a
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,]))
		int_qua_a[i,j]=I
	}
}

#Marginalization of posterior_quad over b
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]))
		int_qua_b[i,j]=I
	}
}
#Marginalization of posterior_quad over m
int_qua_m=array(0, c(dim_a, dim_b))
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
	}
}

#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,]))
		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]))
		int_lin_m[i]=I
}

#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

# 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))
print(strqua)
text(0,2.8,strqua)
text(0,2.1,strlin)
text(0,1.6,strcon)

SS2012-StatisticalMethods.git

new 2.3