#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_con),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 (spinefun) is computationally expensive compared to the simple numerical integration by summing up. 


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. Here i'm not sure how to plot the error bars so I plot them like done below. An explanation about how to plot error bar in the tutorial would be nice.
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.