Code snippets from Introductory Stats Lecture, Introduction Likelihood, comments in German

Teaching-Stats-Likelihood-Intro-GER

# Praktische Vorlesung, Einführung Statistik, WS 13/14
# Florian Hartig, http://www.biom.uni-freiburg.de/mitarbeiter/hartig

?read.csv

# Wiederholung likelihood

# Als estes erstellen wir einen leeren Plot

plot(NULL, NULL, xlim=c(-4,6), ylim = c(0,0.5), ylab = "Wahrscheinlichkeit(sdichte)", xlab = "Observed value")

# OK, jetzt nehmen wir mal an ich beobachte nur einen Punkt, mit einem Wert von 1, ich plotte das

points(1,0, cex = 3, pch = 4, col = "red", lwd = 4)

# nehmen wir mal an ich weiß dass diese Beobachtung aus einem Prozess stammt der einen festen Mittelwert hat,
# und wenn man diesen Prozess beobachtet streuen die beobachteten Werte normalverteilt um den Mittelwert mit einer 
# Standardabweichung von 0.5. Meine Frage ist: was ist der wahrscheinlichste Wert für den Mittelwert dieser Normalverteilung?

# zur illustration plotte ich zwei alternative Normalverteilungen

# ein mal mit Mittelwert 0
# hierfür nutze ich die funktion "curve" die direkt eine andere Funktion plottet
curve(dnorm(x, mean = 0, sd = 1 ), add = T, lty = 2, col = "darkred")

# und einmal mit Mittelwert 1.5
# zur Demonstration diesmal ohne curve

xWerte <- seq(-4,6, length.out = 100)
yWerte <- dnorm(xWerte, mean = 1.5, sd = 1 )
lines(xWerte, yWerte, lty =3, col = "darkgreen")

# nebenbei, wie man plotsymbole, farben und Linientyp ändert ist auch noch mal hier http://www.statmethods.net/advgraphs/parameters.html erklärt


# ok, die Likelihood ist definiert als die Wahrscheinlichkeitsdichte der Modelle die wir anschauen an dem beobachteten Punkt


abline(v = 1)

# also, die zweite Funktion hat eine höhere Likelihood, bzw. mean = 1.5 hat eine höhere Likelihood als mean = 0

# die Maximum Likelihood Methode ist einfach alle Parameterwerte durchzugehen und den mit der höchsten Likelihood zu nehmen

# wir können uns das plotten

curve(dnorm(1, mean = x, sd = 1 ), from = -2, to = 4, ylab = "likelihood", xlab = "Parameterwert für den mean")


# natürlich müssen wir das nicht immer mit hand machen, R hat eine Funktion die uns den besten wert gibt

library(MASS)
fitdistr(1, "normal", sd = 1)


# zur Übung, und für kompliziertere Probleme können wir das auch mit Hand machen

# Diese funktion gibt uns den Likelihood Wert für einen gegebenen Parameter
likelihood<- function(x){
  return(-dnorm(1,x,sd=1, log = T))
} 

# die optiom funktion sucht das maximum einer gegeben Funktion
optim(par = 1.5, fn = likelihood, method = "Brent", lower = -1, upper = 3)


# Achtung, in der funktion Likelihood geben wir die negative log likelihood zurück
# das ist üblich, weil man damit numerisch besser arbeiten kann 

# was ist wenn wir mehrere Datenpunkte haben?
# wenn die Daten punkte unabhängig, aber aus der gleichen Verteilung sind
# (dass nennt man übrigens iid = independent, identically distributed)
# dann ist alles ganz einfach ... die Wahrscheinlichkeine multiplizieren sich

# L = L1 * L2

# und wenn man mit den Logarithment arbeitet, dann muss man sein Schulwissen über log rauskramen

# Log(A*B) = Log(A) + Log(B) 
# Log(A/B) = Log(A) - Log(B)

# Also, wenn man mit logarithmierten Likleihoods arbeitet wird alles einfach addiert

# nLL = nLL1 + nLL2

# ok, dann probieren wir das mal, stellen wir uns vor wir ziehen 5x aus einer Normlaverteilung mit Mittelwert 1.34



data <- rnorm(30,mean=1.34)

# alternative data
#data <- rgamma(100,3,1)
#hist(data)

# was sagt eigentlich dieser p-value --> die Wahrscheinlichkeite diese oder extremere Daten zu sehen, gegeben die Null Hypothese
# klar, hier ist die Null Hypothese dass die Verteilung normal ist, aber welche Normalverteilung genau
# --> die verschiedenen Tests unterscheiden sich in der exakten Formulierung der Null-Hypothese, deshalb gibt es auch
# unterschiedliche p-Werte

# wenn man die Daten  signifikant von der Normalverteilung abweichen, d.h. p>0.05, was heist das, sind wir uns sicher dass der
# zu Grunde liegende Prozess nicht normalverteilt ist?

# Antwort: bei Signifikanzlevel alpha = 0.05 erwarten wir in 5% "falsche Positive" (Typ I Fehler)



plot(NULL, NULL, xlim=c(-2,5), ylim = c(0,0.5), ylab = "Wahrscheinlichkeit(sdichte)", xlab = "Observed value")

points(data,rep(0,length(data)), cex = 3, pch = 3, col = "red", lwd = 1)

result <- fitdistr(data, "normal")
result


# estimate
curve(dnorm(x, mean = result$estimate[1], sd = result$estimate[2] ), add = T, lty = 2, col = "darkred")

# true distribution
curve(dnorm(x, mean = 1.34, sd = 1 ), add = T, lty = 4, col = "darkblue")


# das ganze auch noch mal per Hand

# Diese funktion gibt uns den Likelihood Wert für einen gegebenen Parameter
likelihood<- function(x){
  
  singleLikelihoods = -dnorm(data,x[1],sd=x[2], log = T)
  
  totalLikelihood = sum(singleLikelihoods)
  
  return(totalLikelihood)
} 

# die optiom funktion sucht das maximum einer gegeben Funktion
optim(par = c(1.5,1.2), fn = likelihood)




# vergleiche von verschiedenen Funktionen 



# Testen ob die Daten auch wirklich normal sind

library(nortest)
lillie.test(data)

shapiro.test(data)




# graphische methoden

plot(ecdf(data), pch="", verticals=T, las=1, main="Normalverteilung vs Daten")

curve(pnorm(x, mean=result$estimate[1], sd=result$estimate[2]), add=T, lwd=3, col="red")


# gebräuchlicher: qq plots um die erwarteten und beobachteten Quantile zu vergleichen
qqnorm(data)
qqline(data, col = "red")

# das geht auch mit anderen funktionen über ?ppoints




# finally, comparison via LL and AIC
# remember AIC = -2*log-likelihood + k*npar, where npar represents the number of parameters in the fitted model, and k = 2 for the usual AIC, or k = log(n) (n being the number of observations) for the so-called BIC or SBC (Schwarz's Bayesian criterion).

logLik(result)

AIC(result)


result2 <- fitdistr(data, "gamma")
result2

logLik(result2)
AIC(result2)