#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) {
	#summand = (data[i,2] - m(data[i,1], Omegam, M))**2/data[i,3]**2
	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)
#Erzeugen der zwei Vektoren für schrittweises Wachsen von Omegam und M
Omegam_vec=seq(0.0,1.0, by=0.05)
M_vec=seq(15.5, 16.6, 0.01)
Omegam_vec
M_vec

#Startwerte
chi_min=999999999
Omegam_min=0
M_min=0

#zweifach verschachtelte Schleife für Omegam und M
chisquare_min = chisquare(0,0) #set initial value
for (i in (1: length(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)


jmax=length(Omegam_vec)
kmax=length(M_vec)
j=1
k=1
# generate a posterior matrix for the contour plot
postplot = matrix(nrow=jmax,ncol=kmax)
while (j<=jmax){ k=1
	while (k<=kmax) {
	postplot[j,k]=posterior(Omegam_vec[j],M_vec[k])
	k=k+1}
	j=j+1}
quartz()
contour(Omegam_vec,M_vec,postplot,drawlabels=FALSE,xlab='Omega',ylab='M')


Mmarg=function(Ai) {I=integrate(posterior,13,18,A=Ai)
	return(I$value)
	}