#---------------# # Loesung 7 # #---------------# # Aufg 1b) N = 1000000 x = rnorm(N) y = rnorm(N) z = rnorm(N) F1 = 2*pi*(x^2+y^2) F2 = F1*ifelse(x^2+y^2<1,1,0) F3 = (2*pi)^(3/2) * sqrt(x^2+y^2+z^2)*ifelse(x^2+y^2+z^2<1,1,0) I1 = sum(F1)/N I2 = sum(F2)/N I3 = sum(F3)/N I1 4*pi I2 4*pi*(1-3/2*exp(-1/2)) I3 8*pi*(1-3/2*exp(-1/2)) # Aufg 1c) N = 1000000 x = runif(N,-1,1) y = runif(N,-1,1) z = runif(N,-1,1) F2 = 2*2*(x^2+y^2) * exp(-(x^2+y^2)/2) * ifelse( x^2+y^2<1 , 1 , 0 ) F3 = 2*2*2*sqrt(x^2+y^2+z^2) * exp(-(x^2+y^2+z^2)/2) * ifelse( x^2+y^2+z^2<1 , 1 , 0 ) res2 = sum(F2)/N res3 = sum(F3)/N res2 4*pi*(1-3/2*exp(-1/2)) res3 8*pi*(1-3/2*exp(-1/2)) # Aufg 2 n = 5 N = 100 L = rep(0,20) ell1 = rep(0,20) ell2 = rep(0,20) for(i in 1:30) # wir machen 30 Mal eine Maximum-Likelihood-Schaetzung { x = rt(N,n) for(k in 1:20) { L[k] = prod(dt(x,k)) ell1[k] = log(prod(dt(x,k))) ell2[k] = sum(log(dt(x,k))) } n1 = which.max(L) n2 = which.max(ell1) n3 = which.max(ell2) print(c(n1,n2,n3)) par(mfrow=c(1,2)) plot(L) plot(ell1) points(ell2,col="red") readline("press enter to see next realization..") } # Jetzt mit N=1000 Zufallszahlen: n = 5 N = 1000 L = rep(0,20) ell1 = rep(0,20) ell2 = rep(0,20) for(i in 1:30) { x = rt(N,n) for(k in 1:20) { L[k] = prod(dt(x,k)) # das ist theoretisch ungleich 0, # aber numerisch eine 0 ell1[k] = log(prod(dt(x,k))) # das ist dann also log(0) = -Infinity ell2[k] = sum(log(dt(x,k))) # das ist okay } n1 = which.max(L) n2 = which.max(ell1) n3 = which.max(ell2) print(c(n1,n2,n3)) par(mfrow=c(1,2)) plot(L) plot(ell2) readline("press enter to see next realization..") } ell1 # -Inf