#Problem 1
library(Hmisc)
#definition of a straight line function
linrel=function(x,m,b) m*x+b
#sigvar = 0.4
sigvar = 1
#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 =  sigvar)
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/(sigvar)**2)
chi2_lin = function(m,b) sum((mock_vec - linrel(mock_vec,m,b))**2/(sigvar)**2)
chi2_con = function(b) sum((mock_vec - conrel(b))**2/(sigvar)**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 and m
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, (posterior_qua)[j,i,])
		int_qua_a[i,j]=I
	}
}

int_qua_am=array(0, c(dim_b))
for (i in (1:dim_b)){
		I=int(m_vec, (int_qua_a)[,i])
		int_qua_am[i]=I
}
#Marginalization of posterior_quad over b and m
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, (posterior_qua)[,i,j])
		int_qua_b[i,j]=I
	}
}

int_qua_bm=array(0, c(dim_a))
for (i in (1:dim_a)){
		I=int(m_vec, (int_qua_b)[,i])
		int_qua_bm[i]=I
}

#Marginalization of posterior_quad over b and a

int_qua_ba=array(0, c(dim_m))
for (i in (1:dim_m)){
		I=int(b_vec, (int_qua_a)[,i])
		int_qua_ba[i]=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, (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, (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, (posterior_con))
int_con_b=posterior_con

#search for best sigma
sdtest=seq(0.01,1,by=0.01)

sigma_test=function(A,mu,sd,xvec,f){
	chi=0
	for(j in (1:length(xvec)))
		{
			chi=chi+(f[j]-(A*exp(-0.5*(((xvec[j]-mu)**2)/(sd**2))))/sigvar)**2
        }
        return(chi)
    }

# create grids for vaious sigma functions
# later, find min of grid

grid_con_b <-array(0, dim=c(length(sdtest)))
grid_lin_b <-array(0, dim=c(length(sdtest)))
grid_lin_m <-array(0, dim=c(length(sdtest)))
grid_qua_b <-array(0, dim=c(length(sdtest)))
grid_qua_m <-array(0, dim=c(length(sdtest)))
grid_qua_a <-array(0, dim=c(length(sdtest)))

for(i in (1:length(sdtest)))
{
	grid_con_b[i] = sigma_test(max(int_con_b),best_b_con,sdtest[i], b_vec,int_con_b)
	grid_lin_b[i] = sigma_test(max(int_lin_b),best_b_lin,sdtest[i], b_vec,int_lin_b)
	grid_lin_m[i] = sigma_test(max(int_lin_m),best_m_lin,sdtest[i], m_vec,int_lin_m)
	grid_qua_b[i] = sigma_test(max(int_qua_am),best_b_qua,sdtest[i], b_vec,int_qua_am)
	grid_qua_m[i] = sigma_test(max(int_qua_ba),best_m_qua,sdtest[i], m_vec,int_qua_ba)
	grid_qua_a[i] = sigma_test(max(int_qua_bm),best_a_qua,sdtest[i], a_vec,int_qua_bm)
}
err_con_b=sdtest[which(grid_con_b==min(grid_con_b),arr.ind=TRUE)]
err_lin_b=sdtest[which(grid_lin_b==min(grid_lin_b),arr.ind=TRUE)]
err_lin_m=sdtest[which(grid_lin_m==min(grid_lin_m),arr.ind=TRUE)]
err_qua_b=sdtest[which(grid_qua_b==min(grid_qua_b),arr.ind=TRUE)]
err_qua_m=sdtest[which(grid_qua_m==min(grid_qua_m),arr.ind=TRUE)]
err_qua_a=sdtest[which(grid_qua_a==min(grid_qua_a),arr.ind=TRUE)]
err_lin_m = err_lin_m[1]




# The integration via fit is computationally expensive compared to the simple numerical integration 


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))
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))
strcon = paste("constant: b with b=",round(best_b_con,2),"+-",round(err_con_b,2))
print(strqua)
text(0,2.8,strqua)
text(0,2.1,strlin)
text(0,1.6,strcon)

#Plot errorbars quad
for( x in (x_vec_2)){

tmperrup = quarel(x,best_a_qua,best_m_qua,best_b_qua) + (err_qua_a*x**2+x*err_qua_m+err_qua_b)
tmperrdo = quarel(x,best_a_qua,best_m_qua,best_b_qua) - (err_qua_a*x**2+x*err_qua_m+err_qua_b)
errbar(x,quarel(x,best_a_qua,best_m_qua,best_b_qua),tmperrup,tmperrdo,add=TRUE)

tmperrup = linrel(x,best_m_lin,best_b_lin) + (x*err_lin_m+err_lin_b)
tmperrdo = linrel(x,best_m_lin,best_b_lin) - (x*err_lin_m+err_lin_b)
errbar(x,linrel(x,best_m_lin,best_b_lin),tmperrup,tmperrdo,add=TRUE)

tmperrup = best_b_lin + err_con_b
tmperrdo = best_b_lin - err_con_b
#print(tmperrup)
#print(tmperrdo)
#print(x)
#print(best_b_con)
errbar(x,best_b_con,tmperrup,tmperrdo,add=TRUE)
}
#Problem 4
# calculate the evidence
evidence_qua=err_qua_a*err_qua_m*err_qua_b*max(posterior_qua)*(1/sum(posterior_qua))
evidence_lin=err_lin_m*err_lin_b*max(posterior_lin)*(1/sum(posterior_lin))
evidence_con=err_con_b*max(posterior_con)* (1/sum(posterior_con))
print('evidence')
print('quadratic')
print(evidence_qua)
print('linear')
print(evidence_lin)
print('const.')
print(evidence_con)

#the constant would be the best model. This is in contrast to the naive expectation that straight line model would fit best.

#Problem 5
#To run the script with sigma 1 change the variable sigvar to 1
#Theoretically the evidence should point to the quadratic model because it could better fit the new mock data with higher sigvar. We get the constant again as the best model. Forthermor we gat bigger error bars and so a better evidence.