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