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