Statistik II

Maximum-Likelihood Schätzung, AIC/BIC

Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.

Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.

Versionshistory:

• v1.0: erste online-gestellte Version (14.2.2024)
• v1.1: leichte Änderungen bei Einführung von Likelihood und Likelihood Funktion (3.5.2024)

6 Maximum-Likelihood Schätzung, AIC/BIC

In Statistik I hatten wir uns bereits ausgiebig mit Schätzern und deren Gütekriterien beschäftigt. In den bisherigen Fällen haben wir die jeweiligen Schätzer analytisch hergeleitet (mit Ausnahme der Koeffizienten im Fall der multiplen Regression, was wir in Teil 13 von Statistik II nachholen werden). Die wichtigsten Kriterien die wir dabei betrachtet haben waren, dass ein Schätzer (mindestens) erwartungstreu und konsistent sein soll. Es gibt aber noch eine andere Heransgehensweise um Schätzer zu bestimmen, die insbesondere auch dann zum Tragen kommt, wenn es keine analytisch bestimmbare Lösung gibt: die Maximum-Likelihood (ML) Schätzung.

# Pakete die in diesem Teil benutzt werden:
library(MASS)
library(dfoptim)

6.1 Maximum-Likelihood Schätzung von Parametern

Wir sind bereits des Öfteren Fragen begegnet wie “Mit welcher Wahrscheinlichkeit treten Werte \(\geq 100\) auf?”, wenn angenommen wird, die zugrunde liegende Variable sei normalverteilt mit \(\mu = 50\) und \(\sigma = 20\). Hierbei sind die Parameter der Verteilung bereits bekannt und wir fragen nach der Wahrscheinlichkeit von Daten. Allerdings sind in vielen (Forschungs-)Situationen genau diese Parameter unbekannt und sie müssen auf Basis von Daten (d.h. den vorliegenden Werten einer Stichprobe) bestimmt werden. Dafür benötigen wir die ML Methode und die resultierenden Schätzer für diese Parameter heißen dann ML Schätzer.

ML Schätzung kann als eine Methode aufgefasst werden, mit der (Populations-)Parameter verschiedenster Verteilungen so aus Daten geschätzt werden, dass die “Wahrscheinlichkeit” – oder genauer: die Likelihood – der beobachteten Daten maximiert wird.

6.1.1 Beispieldaten

Als Beispiel für die folgenden Argumentationen und Berechnungen verwenden wir die gemessenen Größen (in Zentimetern) von \(n=10\) Personen. Da wir die Daten hier simulieren, wissen wir bereits, dass diese aus einer Normalverteilung mit \(\mu = 170\) und \(\sigma=20\) stammen. Natürlich ist in der generierten Stichprobe der Mittelwert \(M\) und die Standardabweichung \(S\) (bzw. \(\hat{S}\)) von diesen Parametern aber abweichend:

set.seed(1)
data <- round(rnorm(n = 20,
                    mean = 170,
                    sd = 20))
data

##  [1] 157 174 153 202 177 154 180 185 182 164 200 178 158 126 192 169 170 189 186
## [20] 182

mean(data)               # Mittelwert

## [1] 173.9

n <- length(data)        # Anzahl Datenpunkte... für:
sqrt(var(data)*(n-1)/n)  # Stichproben-Standardabweichung

## [1] 17.7395

6.1.2 Grundlagen

Im Folgenden bezeichnen wir mit \(X\) die empirischen Daten \(x_1,\ldots,x_n\) der \(n\)-vielen Versuchspersonen. Den Satz an Parametern bezeichnen wir mit \(\theta\) (einem kleinen Theta), im Fall der Normalverteilung also \(\theta = (\mu,\sigma\)).

Das Ziel wird sein, diejenige Dichtefunktion (oder Wahrscheinlichkeitsfunktion) \(f(X;\theta)\) aus allen möglichen Parametrisierungen zu bestimmen, die die beobachteten Daten \(X\) am “wahrscheinlichsten” produziert hat. Auf dem Weg dahin machen wir uns zunächst eine Sache (noch einmal) klar. Die Normalverteilung wird üblicherweise geschrieben als: \[ f(x; \mu, \sigma)=\frac{1}{\sqrt{2\pi\sigma^2}}\cdot e^{-\frac{(x-\mu)^2}{2\sigma^2}} \] Streng genommen ist sie also eine Funktion von drei Variablen, die allesamt variiert werden können und dann den entsprechenden Funktionswert \(f\) liefert. Was wir aber üblicherweise getan haben, ist sie als Funktion von \(x\) zu betrachten während \(\mu\) und \(\sigma\) feststehen. In der folgenden Abbildung ist \(\mu=2\) und \(\sigma=1.5\) und wir fragen: Welchen Wert nimmt die Normalverteilung dann für z.B. \(x=1\) an. Die Antwort ist dann \(f(1)=0.213\):

Nun gibt es den ersten wichtigen Begriff: Den Wert, welchen wir durch die Evaluation der Dichtefunktion für \(x\), \(\mu\) und \(\sigma\) erhalten, bezeichnen wir als “Likelihood” und er gibt an wie “plausibel” ein Datenpunkt ist.

Wichtig ist dabei anzumerken, dass die Likelihood keine Wahrscheinlichkeit ist. Wahrscheinlichkeiten sind (im stetigen Fall) immer noch Flächen unter Dichtefunktionen und die Wahrscheinlichkeit eines konkreten Wertes \(x\) ist bei einer stetigen Verteilung stets 0. Dennoch ist es bspw. in der letzten Abbildung relativ betrachtet “plausibler” einen Wert um 2 zu ziehen als einen Wert um 1. Tatsächlich könnte ein Intervall um einen Wert \(a\) gelegt werden mit \([a\pm\epsilon]\) (für ein beliebig kleines \(\epsilon\)) und die Wahrscheinlichkeit einen Wert in diesem Bereich zu ziehen wäre bei \(a=2\) größer als bei \(a=1\).

6.1.3 Likelihood-Funktion

Nun drehen wir die Perspektive und betrachten unsere Daten als fixiert bzw. gegeben. Die Dichtefunktion \(f\) ist dann nur noch eine Funktion von \(\theta\) (z.B. von \(\mu\) und/oder \(\sigma\)).

Eine solche Funktion nennen wir Likelihood-Funktion: \[ \mathcal{L}(\theta|X) = f(X;\theta) \] Wichtig ist hierbei zu verstehen, dass die Likelihood-Funktion nichts anderes ist als eine Dichtefunktion, unter Annahme, dass \(X\) bekannt wäre! Wir greifen also bei der Likelihood-Funktion auch auf die Gleichung der Dichtefunktion zurück. Damit wird die Bestimmung in R mittels der d-Funktionen möglich, z.B. dnorm() im Fall der Normalverteilung, dt() im Fall der t-Verteilung oder dbinom() im Fall der Binomialverteilung.

Berechnen wir nun einen ersten Likelihood-Wert. Der erste Wert in unserem Beispieldatensatz ist \(x_1=157\) und wir nehmen an (bzw. wissen in unserem Fall), dass die Normalverteilung die korrekte zugrunde liegende Verteilung ist. Da wir die Verteilung der Normalverteilung ausschreiben können, können wir also \(\mathcal{L}(157)\), zunächst mit den Parameter \(\mu = 170\) und \(\sigma=20\) (also denen, die wir zum Produzieren der Daten verwendet haben), wie folgt berechnen: \[ \begin{aligned} f(157; \mu = 170, \sigma=20)&=\frac{1}{\sqrt{2\pi\sigma^2}}\cdot e^{-\frac{(x-\mu)^2}{2\sigma^2}}\\ f(157; \mu = 170, \sigma=20)&=\frac{1}{\sqrt{2\pi 20^2}}\cdot e^{-\frac{(157-170)^2}{2\cdot 20^2}}\\ &=0.01994711\cdot 2.718282^{-\frac{169}{800}}\\ &=0.01994711\cdot 2.718282^{-0.21125}\\ &=0.01614861 \end{aligned} \]

Den gleichen Wert (abgesehen von Rundungsungenauigkeiten in der obigen Rechnung) erhalten wir, wenn wir die Funktion dnorm() entsprechend verwenden:

dnorm(x = 157,
      mean = 170,
      sd = 20)

## [1] 0.01614862

Verwenden wir andere Parameter für \(\theta=(\mu,\sigma)\), dann erhalten wir entsprechend auch einen anderen (hier kleineren) Likelihood-Wert. Der Wert stammt also “eher” aus der gerade verwendeten Verteilung mit \(\theta=(170,20)\) als aus der im folgenden Beispiel mit \(\theta=(100,20)\) (was hier trivial ist, da wir ja gerade die “korrekten” Verteilungsparameter verwendet haben):

dnorm(x = 157,
      mean = 100,
      sd = 20)

## [1] 0.0003436383

Der nächste Schritt ist nun die Bestimmung der Likelihood-Funktion unter Betrachtung der gesamten Stichprobe gleichzeitig. Jeder einzelne Datenpunkt kann als Ereignis aufgefasst werden und unter der Annahme, alle diese Ereignisse seien unabhängig voneinander, ergibt sich die Likelihood-Funktion als Produkt der einzelnen Likelihoods: \[ \begin{aligned} \mathcal{L}(\theta|X) &= f(X;\theta)\\ &=f(x_1,x_2,\ldots,x_n;\theta)\\ &=f(x_1;\theta)\cdot f(x_2;\theta)\cdot \ldots f(x_n;\theta)\\ &=\prod_{i=1}^n f(x_i;\theta) \end{aligned} \] Das Maximum dieser Funktion in Abhängigkeit von \(\theta\) kann nun gefunden werden, indem die (partielle) Ableitung bzgl. \(\theta\) gleich Null gesetzt wird. Technisch gibt es hier aber noch zwei Herausforderungen: Zum Einen erfordern Ableitungen von Funktionen die aus Produkten bestehen die etwas unhandliche Kettenregel. Zum Anderen sind die einzelnen Likelihoods oft (aber nicht immer) kleiner als Eins und die Gesamt-Likelihood \(\mathcal{L}\) wird dadurch extrem klein. Dies bereitet Computern tatsächlich ein Problem, da diese nur eine bestimmte Anzahl an Dezimalstellen repräsentieren können. Üblicher ist daher die sogenannte logarithmierte Likelihood oder log-Likelihood Funktion zu verwenden.

6.1.4 Logarithmierte Likelihood Funktion (log-Likelihood)

Durch die Logarithmierung werden zunächst “kleine” Zahlen wieder betragsmäßig größer:

log(0.0001)

## [1] -9.21034

log(0.000000001)

## [1] -20.72327

Darüber hinaus erlaubt die Produktregel beim Rechnen mit Logarithmen uns, das Produkt in eine Summe umzuformulieren, da sie besagt: \[\log_b(y\cdot z) = \log_b(y)+\log_b(z)\]

# ein Beispiel dazu
y <- 2
z <- 3
log(y*z)

## [1] 1.791759

log(y) + log(z)

## [1] 1.791759

Wenngleich die Basis des verwendeten Logarithmus eigentlich egal ist, wird in der Regel der natürliche Logarithmus zur Basis \(e\) verwendet. Umgeschrieben ergibt sich die log-Likelihood \(\ell\) daher als: \[ \begin{aligned} \ell = \ln(\mathcal{L}) &= \ln\left[ \prod_{i=1}^n f(x_i;\theta) \right] \\ &= \sum_{i=1}^n \ln \left[ f(x_i;\theta) \right]\\ &= \ln[f(x_1;\theta)]+\ln[f(x_2;\theta)]+\ldots+\ln[f(x_n;\theta)] \end{aligned} \]

Das Schöne an dieser Transformation ist, dass das Maximum von Likelihood und log-Likehood gleich ist. Die folgende Abbildung stellt links die Likelihood und rechts die log-Likelihood für verschiedene Werte von \(\mu\) (bei konstantem \(\sigma=20\)) anhand der oben generierten Daten dar (man beachte hierbei die Skalierungen der \(y\)-Achsen). Zusätzlich ist als senkrechte Linie das jeweilige Maximum eingezeichnet.

Anmerkung: Der Wert der log-Likelihood Funktion wird umgangssprachlich oft auch nur log-Likelihood genannt.

6.2 Bestimmung von Maximum-Likelihood Schätzern

Nun da wir die log-Likelihood Funktion kennen, stellt sich die Frage: Wie finden wir deren Maximum praktisch? Hier gibt es zwei Antworten. In manchen Fällen, lassen sich analytische Lösungen finden, die dann als Gleichung geschrieben werden können. In (vielen) anderen Fällen hingegen ist dies nicht möglich und wir müssen auf numerische Verfahren zurückgreifen, die verschiedenste Parameterkonstellationen durchgehen und das Maximum der log-Likelihood so bestimmen.

6.2.1 Analytische Lösung

Analytische Lösungen betrachten im Wesentlichen die (partiellen) Ableitungen der log-Likelihood Funktion bzgl. der Parameter in \(\theta\) und setzen diese gleich Null (Dies wird manchmal als Score-Gleichung oder Score-Funktion bezeichnet). Tatsächlich können die ML Schätzer der Parameter der Normalverteilung analytisch bestimmt werden und wir beginnen mit dem Parameter \(\mu\).

6.2.1.1 Der Parameter \(\mu\) der Normalverteilung

Dazu beginnen wir mit der Likelihood \(\mathcal{L}\) und formen diese unter Verwendung verschiedener Rechenregeln etwas um. Zunächst gilt \(e^y\cdot e^z=e^{y+z}\) bzw. in Verallgemeinerung: \[ \prod_{i=1}^n e^{x_i}=e^{\sum_{i=1}^n x_i} \] Darüber hinaus gilt für eine Konstante \(z\) \[ \prod_{i=1}^nz=z^n \] und für Produkte zweier Variablen gilt \[ \prod_{i=1}^nx_iy_i=\prod_{i=1}^nx_i\cdot \prod_{i=1}^ny_i \] Diese Regeln werden im Folgenden angewendet, um die Likelihood-Funktion etwas umzuschreiben: \[ \begin{aligned} \mathcal{L}&=\prod_{i=1}^n\left[\frac{1}{\sqrt{2\pi\sigma^2}}\cdot e^{-\frac{(x_i-\mu)^2}{2\sigma^2}}\right]\\ &=\prod_{i=1}^n\frac{1}{\sqrt{2\pi\sigma^2}}\cdot \prod_{i=1}^ne^{-\frac{(x_i-\mu)^2}{2\sigma^2}}\\ &=\left[ \frac{1}{\sqrt{2\pi\sigma^2}}\right]^n\cdot \prod_{i=1}^ne^{-\frac{(x_i-\mu)^2}{2\sigma^2}}\\ &=\frac{1}{(\sigma\sqrt{2\pi})^n} \cdot e^{-\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n(x_i-\mu)^2} \end{aligned} \]

Durch Logarithmierung gelangen wir nun zur log-Likelihood-Funktion \(\ell\): \[ \begin{aligned} \ell = \ln[\mathcal{L}]&=\ln\left[ \frac{1}{(\sigma\sqrt{2\pi})^n} \cdot e^{-\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n(x_i-\mu)^2} \right] \\ &=\ln\left[ \frac{1}{(\sigma\sqrt{2\pi})^n} \right] - {\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n(x_i-\mu)^2} \\ &=\ln\left[\sigma^{-n}(2\pi)^{-\frac{n}{2}} \right] - {\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n(x_i-\mu)^2} \\ &=-n\cdot \ln[\sigma]-\frac{n}{2}\cdot\ln[2\pi] - {\frac{1}{2\sigma^2}\cdot\sum_{i=1}^n(x_i-\mu)^2} \\ \end{aligned} \] Nun leiten wir \(\ell\) bzgl. \(\mu\) (partiell) ab und setzen das Ergebnis gleich Null: \[ \frac{\partial \ell}{\partial\mu}=\frac{1}{\sigma^2}\sum_{i=1}^n (x_i-\mu)=0 \] Da \(\sigma>0\) sein muss, kann der Ausdruck nur Null werden, wenn der zweite Teil (also die Summe) gleich Null ist und diese Gleichung kann dann nach \(\mu\) aufgelöst werden: \[ \begin{aligned} &\sum_{i=1}^n (x_i-\mu)=0\\ \Leftrightarrow &\sum_{i=1}^nx_i-\sum_{i=1}^n\mu=0\\ \Leftrightarrow &\sum_{i=1}^nx_i-n\mu=0\\ \Leftrightarrow &\sum_{i=1}^nx_i=n\mu\\ \Leftrightarrow &\frac{1}{n}\sum_{i=1}^nx_i=\mu\\ \end{aligned} \] Diese letzte Formel sollte nun vertraut vorkommen, sie ist nämlich nichts anderes als der Mittelwert, den wir auch in Statistik I schon als erwartungstreuen und konsistenten Schätzer für den Erwartungswert \(\mu\) kennengelernt haben. Das heißt, der Mittelwert ist (auch) der ML-Schätzer für den (Verteilungs-)Parameter \(\mu\) der Normalverteilung.

6.2.1.2 Der Parameter \(\sigma\) der Normalverteilung

Zur Bestimmung des ML-Schätzers für \(\sigma\) leiten wir \(\ell\) entsprechend (partiell) nach \(\sigma\) ab und setzen das Ergebnis gleich Null:

\[ \frac{\partial \ell}{\partial\sigma}=\frac{-n}{\sigma}+\frac{\sum_{i=1}^n (x_i-\mu)^2}{\sigma^3}=0 \] Nun multiplizieren wir beide Seiten mit \(\sigma^3\) und lösen dann nach \(\sigma^2\) auf: \[ \begin{aligned} &\frac{-n}{\sigma}+\frac{\sum_{i=1}^n (x_i-\mu)^2}{\sigma^3}=0\\ \Leftrightarrow &-n\sigma^2+\sum_{i=1}^n (x_i-\mu)^2=0\\ \Leftrightarrow &\sum_{i=1}^n (x_i-\mu)^2=n\sigma^2\\ \Leftrightarrow &\frac{1}{n}\sum_{i=1}^n (x_i-\mu)^2=\sigma^2\\ \end{aligned} \] Auch diese Formel sollte bekannt vorkommen: Sie ist die Formel der Stichprobenvarianz und die Wurzel aus dem Ergebnis die Standardabweichung. Das heißt, die Standardabweichung ist der ML-Schätzer für den (Verteilungs-)Parameter \(\sigma\) der Normalverteilung.

6.2.2 Numerische Lösung

Gerade haben wir analytische Lösungen aufgezeigt, um an ML-Schätzer für \(\mu\) und \(\sigma\) zu kommen. Wie bereits erwähnt ist das häufig nicht möglich. Dann werden numerische Verfahren verwendet, die systematisch nach denjenigen Parameterwerten suchen, die eine Funktion maximieren. Typische Funktionen, die dies in R tun können sind optim() und nmkb() (aus dem Paket dfoptim).

Derartige Funktionen benötigen zunächst eine “Kostenfunktion”, die denjenigen Wert berechnet der maximiert werden soll. Dazu schreiben wir eine R-Funktion, die für einen gegebenen Datenvektor X und einen Parametervektor \(\theta\) die log-Likelihood berechnet, wobei wir hier wieder \(\mu\) und \(\sigma\) der Normalverteilung für die oben generierten Daten bestimmen wollen:

# Kostenfunktion: Berechnung der log-Likelihood
logL <- function(X, theta) {      # theta wird ein Vektor mit zwei Elementen sein
  mu <- theta[1]                  # mu extrahieren
  sigma <- theta[2]               # sigma extrahieren
  logL <- sum(dnorm(x = X,        # empirische Daten
                    mean = mu,    # Parameter mu
                    sd = sigma,   # Parameter sigma
                    log = TRUE))  # liefert direkt ln(f(x_i)) zurück
  return(logL)
}

Die Funktion nmkb() benötigt nun vor allem Startwerte für die Parameter und sie kann den möglichen Wertebereich auch direkt begrenzen, was hier sinnvoll ist, da z.B. Werte \(\sigma\leq0\) nicht in Frage kommen. Dazu legen wir drei Vektoren an:

start.theta <- c(150,30)  # Startwerte für mu und sigma
lower.theta <- c(100,0)   # untere Grenze für die Parameter
upper.theta <- c(200,50)  # obere Grenze für die Parameter

Nun können wir die Funktion selber aufrufen und das Ergebnis speichern wir im Objekt result:

result <- nmkb(par = start.theta,
               lower = lower.theta,
               upper = upper.theta,
               fn = logL,
               X = data,
               control = list(maximize = TRUE))   # da maximiert werden soll

Das Objekt result enthält eine ganze Reihe von Werten, besonders interessant sind für uns nun diejenigen Parameter, die zur größten log-Likelihood geführt haben, die wir direkt mit dem Mittelwert und der Standardabweichung der Daten vergleichen:

result$par               # beste Parameter

## [1] 173.90056  17.74191

mean(data)               # Mittelwert

## [1] 173.9

n <- length(data)        # Anzahl Datenpunkte... für:
sqrt(var(data)*(n-1)/n)  # Stichproben-Standardabweichung

## [1] 17.7395

Wir erhalten durch den Algorithmus also minimal unterschiedliche Ergebnisse im Vergleich zu den analytisch bestimmten Werten. Dies liegt aber daran, dass der Algorithmus das gesuchte Maximum nicht hinreichend genau abtastet bzw. die minimale Schrittbreite (standardmäßig) nicht fein eingestellt genug ist. Verfeinern wir dies durch Änderung der entsprechenden Parameter der Funktion, dann erhalten wir die gleichen oder nur noch minimal abweichende Werte. Kleine Abweichungen bleiben allerdings, da ein Computer eben nur eine begrenzte Zahl von Dezimalstellen repräsentieren kann.

In Abschnitt 4 werden wir noch eine Funktion kennenlernen, die die Suche nach den besten Parametern – das “Fitten” einer Verteilung an Daten – sehr einfach gestaltet.

6.2.3 Vergleich von Schätzern bzw. Schätzmethoden

Die nachfolgende Tabelle stellt die Eigenschaften der betrachteten Schätzer gegenüber. Ersichtlich ist, dass die “besten” Schätzer sich unterscheiden können in Abhängigkeit von der verwendeten Konstruktionsmethode. Allerdings gilt für ML-Schätzer, dass diese asymptotisch erwartungstreu sind, d.h. \(E(T)=\tau\) wenn \(n\rightarrow \infty\).

6.2.4 Konfidenzintervalle für ML-Schätzer

Ein großer Vorteil von ML-Schätzern ist es, dass in fast allen Fällen Konfidenzintervalle berechnet werden können. Hierfür ist die negative zweite (partielle) Ableitung der log-Likelihood Funktion relevant, die auch als Fisher-Information bezeichnet wird. Die Fisher-Information für das Maximum der log-Likelihood Funktion (also den Maximum-Likelihood Schätzer) nennt man beobachtete Fisher-Information.

Bezeichnen wir mit \(\hat{\theta}\) den Maximum-Likelihood Schätzer für \(\theta\) und mit \(I\) die Fisher-Information, dann kann man herleiten, dass für den Standardfehler von \(\hat{\theta}\) gilt: \[ SE_{\hat{\theta}}=\frac{1}{\sqrt{I(\hat{\theta})}} \] Mit anderen Worten ist also die Wurzel aus der inversen beobachteten Fisher-Information ein Standardfehler von \(\hat{\theta}\). Daraus wiederum lassen sich Konfidenzintervalle berechnen.

6.3 Vergleich von Fits

Es gibt viele Situationen, wo Daten prinzipiell durch mehrere Verteilungen beschrieben werden könnten und es stellt sich dann die Frage: Welche Verteilung lässt sich besser an die Daten fitten? Oder anders: Welche Verteilung ist plausibler als Quelle der Daten?

6.3.1 Log-Likelihood

Mit der log-Likelihood \(\ell\) haben wir direkt ein Maß, welches hierfür herangezogen werden kann. Angenommen wir hätten für die obigen Daten auch die sogenannte log-normale Verteilung ausprobiert (die wir aber nicht kennengelernt haben) und erhalten \[\ell_\text{normal}=-100\qquad\text{und}\qquad \ell_\text{log-normal}=-200\] dann bedeutet dies, dass die Normalverteilung besser fittet, also plausibler als Quelle für die Daten ist. Allgemein gilt: je größer die log-Likelihood \(\ell\), desto besser ist der Fit.

6.3.2 Informationstheoretische Werte: AIC und BIC

Ein Problem ergibt sich, da Verteilungen verschieden viele Parameter haben können und eine Verteilung mit mehr Parametern i.d.R. flexibler ist und daher die log-Likelihood i.d.R. größer ist. Daher wurde vorgeschlagen, dass der entsprechende Wert noch für die Anzahl der Parameter korrigiert werden soll. Ein sehr einflussreicher Vorschlag dazu stammt von Akaike (1974) und wird heute Akaike Information Criterion (AIC) genannt: \[AIC=-2\ell+2p\] wobei \(p\) die Anzahl der Parameter darstellt. Für diese Formel gibt es noch verschiedene Korrekturen für die Stichprobengröße, wobei die vielleicht wichtigste Alternative das Bayesian Information Criterion (BIC) ist (Schwarz, 1978): \[BIC = -2\ell+p\log(n)\] wobei (streng genommen) \(n\) der sogenannte “limitierende Stichprobenumfang” ist, was aber bei kontinuierlich verteilten Daten der Anzahl der Datenpunkte \(n\) entspricht.

Da die log-Likelihood mit einem negativen Vorzeichen in die Berechnung beider Maße eingeht, sprechen nun wieder kleinere Werte für einen besseren Fit.

6.3.3 Inferenzstatistische Vergleiche

Auch wenn es eigentlich gegen die Idee des Ansatzes ist, lassen sich auch inferenzstatistische Vergleiche (statt einfacher größer/kleiner-Bewertungen) durchführen. Dabei macht man sich zunutze, dass die sogenannte Deviance-Statistik \[\text{Deviance}=-2\cdot\ell\] approximativ \(\chi^2\)-verteilt ist. Bei genesteten Modellen sind auch Differenzen der jeweiligen Deviance-Werte approximativ \(\chi^2\)-verteilt, was wir im Rahmen der logistischen Regression für Modellvergleiche benötigen werden.

6.4 Praktische Beispiele

Während die bisherigen Abschnitte eher theoretisch ausgelegt waren, werden nun noch zwei Beispiele vorgestellt, in deren Rahmen auch Funktionen verwendet werden, die das Fitten von Verteilungen an Daten automatisch übernehmen, ohne dass wir dazu einen Maximierungsalgorithmus verwenden. Eine Funktion dafür ist fitdistr() aus dem Paket MASS.

6.4.1 Fitten von Daten an Normalverteilung

Wir verwenden abermals die oben generierten Daten data und verwenden nun die Funktion fitdistr() zur Schätzung der Parameter der Normalverteilung:

result <- fitdistr(x = data,             # Datensatz X
                   densfun = "normal")   # welche Verteilung soll gefittet werden?
result$estimate                          # Ausgabe der Parameter

##     mean       sd 
## 173.9000  17.7395

Wie wir sehen, passen die besten Parameter sehr gut zu dem, was wir oben berechnet haben. Tatsächlich stimmen sie in dem Beispiel mit dem Mittelwert und der Standardabweichung überein, da die Funktion intelligent genug ist zu wissen, dass es eine analytische Lösung zur Bestimmung der ML-Schätzer für \(\mu\) und \(\sigma\) gibt…. :)

6.4.2 Vergleich des Fits (z.B. von Reaktionszeiten) zwischen Normal- und Gamma-Verteilung

Nun wollen wir noch ein weiteres Beispiel betrachten. Die folgende Abbildung stellt ein Histogramm empirischer Daten dar, wie sie z.B. von Reaktionszeiten stammen können:

In einem ersten Ansatz könnte man denken, hier würde eine Normalverteilung zugrunde liegen und die entsprechenden Parameter \(\mu\) und \(\sigma\) werden bestimmt:

result.n <- fitdistr(x = data2,            # Datensatz
                     densfun = "normal")   # welche Verteilung soll gefittet werden?
result.n$estimate                          # Ausgabe der Parameter

##      mean        sd 
## 1.0403465 0.7279006

Wir extrahieren die beiden Parameter und zeichnen das Histogramm erneut, allerdings dieses Mal mit der gerade berechneten Normalverteilung hinzugefügt:

mu <- result.n$estimate[1]
sigma <- result.n$estimate[2]
curve(dnorm(x, mean = mu, sd = sigma),  # Dichte der Normalverteilung mit mu und sigma
      from = -2, to = 5,                # von ... bis
      add = TRUE)                       # dem PLot hinzufügen

Hier wird offensichtlich, dass eine Normalverteilung nicht wirklich zu den Daten passt. Als eine Alternative bietet sich z.B. die Gamma-Verteilung an, deren zwei Parameter rate und shape genannt werden. Das Fitten der Gamma-Verteilung an die Daten geht wieder mit der Funktion fitdistr():

result.g <- fitdistr(x = data2,             # Datensatz
                     densfun = "gamma")     # welche Verteilung soll gefittet werden?

## Warning in densfun(x, parm[1], parm[2], ...): NaNs wurden erzeugt

result.g$estimate                           # Ausgabe der Parameter

##    shape     rate 
## 2.028006 1.949357

Auch hier extrahieren wir die Parameter und plotten das Histogramm mit der darübergelegten Gamma-Verteilung:

shape <- result.g$estimate[1]
rate <- result.g$estimate[2]
curve(dgamma(x, shape = shape, rate = rate),  # Dichte der Gamma-Verteilung...
      from = -2, to = 5,
      add = TRUE)

##    shape     rate 
## 2.028006 1.949357

Offenbar passt diese Verteilung optisch deutlich besser zu den Daten. Dies sollte sich natürlich auch in den jeweiligen Werten für die log-Likelihood, \(AIC\) und \(BIC\) widerspiegeln. Zur Berechnung gibt es die Funktionen logLik(), AIC() und BIC, denen jeweils das Ergebnis der Funktion fitdistr() übergeben wird:

# log-Likelihood
logLik(result.n)   # Normalverteilung

## 'log Lik.' -1101.348 (df=2)

logLik(result.g)   # Gamma-Verteilung

## 'log Lik.' -919.5675 (df=2)

# AIC
AIC(result.n)

## [1] 2206.695

AIC(result.g)

## [1] 1843.135

# BIC
BIC(result.n)

## [1] 2216.511

BIC(result.g)

## [1] 1852.95

Auf Basis aller Werte ergibt sich, dass die Gamma-Verteilung besser die Daten fittet als die Normalverteilung es tut.

6.5 Literatur

Schwarz, G. (1978). Estimating the dimension of a model. The annals of statistics, 6(2), 461-464.

Akaike, H. (1974). A new look at the statistical model identification. IEEE transactions on automatic control, 19(6), 716-723.