```#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)')

```