#Exercise 1
#a)
#setwd("C:/Users/Studium/Desktop/Statistische Methoden/Sheet 4")
input=scan(file = "sn_data_riess.dat", what = list(character(), double(), double(), double()), skip=1, multi.line=FALSE)
data=cbind(input[[2]], input[[3]], input[[4]])
data

#b)
#Constants
H0=72.0 #km/s/Mpc
c = 3*10^5 #km/s


#Hubble's law
Hubble=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))


#H^(-1)
Hinv=function(zint,Omegam) 1.0/Hubble(zint,Omegam)


#Luminosity distance
dL=function(z,Omegami){
	dLsol=z
	for (i in (1:length(z))){
		zarg=z[i]
		I=integrate(Hinv,0.0,zarg,Omegam=Omegami)
		dLsol[i] = c*(1+zarg)*I$value
		}
	return(dLsol)
	}

m=function(z, Omegam, M) M + 5*log10(H0 * dL(z, Omegam))

mag_th=NULL
mag_th[186]=1
for (i in (1:186)){
	mag_th[i] = m(data[i,1], 1, data[i,2]) #Omegam = 1 for a flat universe 
	}

#c)
chisquare=function(Omegam, M) {
	chisqu=0
	for (i in (1:186)){
		chisqu=chisqu+(data[i,2] - m(data[i,1], Omegam, M))**2/data[i,3]**2
	}
	return(chisqu)
}

	
#d)
#create vectors Omegam und M
Omegam_vec=seq(0.0,1.0, by=0.05)
M_vec=seq(15.5, 16.6, 0.01)

chisquare_min = chisquare(0,0) #set initial value
Omegam_min=0
M_min=0

#zweifach verschachtelte Schleife f�r Omegam und M

for (i in M_vec){
	for (j in Omegam_vec){
		if (chisquare(j,i) < chisquare_min){ #compare current value with known smallest value 
			chisquare_min = chisquare(j,i)# set new smallest value
			Omegam_min = j 
			M_min = i
		}
	}
}
chisquare_min
Omegam_min
M_min

#e)
posterior = function(Omegam, M){
	p = exp(-0.5*(chisquare(Omegam, M)- chisquare_min))
	return(p)
}

#f)
Atest=seq(0,1,by = 0.05)	
Btest=seq(15.5,16.5,by = 0.01)

jmax=length(Atest)
kmax=length(Btest)
j=1
k=1

postplot = matrix(nrow=jmax,ncol=kmax)
while (j<=jmax){ k=1
	while (k<=kmax) {
	postplot[j,k]=posterior(Atest[j],Btest[k])
	k=k+1}
	j=j+1}
contour(Atest,Btest,postplot,levels=c(1.0,0.8,0.6,0.4,0.2,0.0),drawlabels=FALSE,xlab='Omegam',ylab='M',xlim=c(0,3),ylim=c(0,3))

#g)
pMmarg=function(Mi) {I=integrate(posterior,13.0,18.0,A=Mi)
	return(I$value)
	}

#h)
pMplot=function(Mi) {dumplot=Mi
	for (i in (1:length(Mi))) {
	     Marg=Mi[i]*1.0
	     dumplot[i]=pMmarg(Marg)
	            }
	return(dumplot)}
	
Mplot=seq(0.1,3.0,by = 0.02)

normM=1.0/max(pMplot(Mplot))
plot(Mplot,normM*pMplot(Mplot),type='l',lwd=2,xlab='M',ylab='prob(M|{N_k},I)')