#----------------------------#
#    Zu week9:  Grundidee    #
#     Maximum-Likelihood     #
#----------------------------#


# Wir erzeugen n Zufallszahlen mit Mittelwert mu 
# und Standardabweichung sigma: 

mu = 15 
sigma = 5

n = 1000

set.seed(123456)
z = rnorm(n,mean=mu,sd=sigma)
z

plot(z)
hist(z,breaks=30,prob=TRUE)
curve(dnorm(x,mu,sigma),add=TRUE,col="red")


# Jetzt: Nehmen wir an, wir haben nur den Vektor z gegeben und wissen, 
# dass es normalverteilte Zufallszahlen sind, aber wir kennen die Werte 
# von mu und sigma nicht. 

# Grundlegendes Problem: Wie kann man die Modellparameter, hier mu und sigma, 
# aus den Daten, hier die z1,...,z1000, zurueckgewinnen? 

# Wir muessen eine Likelihood-Funktion aufstellen und maximieren: 

# exakt:

muML = sum(z)/n
muML
sigmaML = sqrt(  sum( (z-muML)^2 )/n  )
sigmaML


# Numerisches Maximieren der Likelihood-Funktion:

# wir waehlen test-mu's und test-sigma's:

mu = seq(from=10,to=20,by=0.05)
nmu = length(mu)
sigma = seq(from=2,to=8,by=0.02)
nsigma = length(sigma)

logL = matrix(0,nrow=nmu,ncol=nsigma)

for(i in 1:nmu)
{
  for(j in 1:nsigma)
  {
    logL[i,j] = -n*log(sigma[j]) - sum( (z-mu[i])^2 )/(2*sigma[j]^2) 
  }
}


contour(mu,sigma,logL, zlim=c(-2500,0), nlevels=400 )


maxlogL =max(logL)
maxlogL

which(logL==maxlogL, arr.ind=TRUE)

mu[102]
sigma[149]