2 Verteilungen in R

Eine große Anzahl an Verteilungen ist bereits mit der Basisverion von R installiert. Weitere Verteilungen können über Pakete, wie z.B. das fGarch-Paket mit installiert werden. Zu jeder Verteilung existiert in R eine Funktion:

  • zum ziehen von Zufallszahlen
  • zur Berechnung der Wahrscheinlichkeitsmasse bzw. -dichte
  • zur Berechnung der kumulativen Verteilungsfunktion
  • zur Bestimmung des \(\alpha\)-Quantils der Verteilung

Die Syntax ist hier bei allen Verteilungen gleich r-, d-, p- und q-Verteilungsname. Am Beispiel der Normalverteilung:

##  [1] -0.0002771885  0.6151002700  1.7785104858  0.7200973972  1.1742840901
##  [6]  0.6049259361 -0.3414518510 -1.7743312435  1.3067025492 -0.6140824878
##  [1] 0.39894227 0.33018152 0.08204493 0.30782967 0.20020567 0.33223716
##  [7] 0.37635094 0.08265630 0.16987807 0.33038812
##  [1] 0.49988942 0.73075575 0.96233997 0.76426749 0.87985937 0.72738590
##  [7] 0.36638172 0.03800419 0.90434312 0.26958041
## [1] -1.6448536 -0.6744898  0.0000000  0.6744898  1.6448536

2.1 Parametrische Verteilungen

Wie der Name impliziert zeichnen sich parametrische Verteilungen dadurch aus, dass Parameter spezifiziert sein müssen, um die Verteilungen zu definieren. In R werden die Parameter innerhalb der Funktionsaufrufe festgelegt. Um bei einer Verteilungsfunktion heraus zu finden, welche Parameter spezifiziert werden müssen bietet es sich sich an in der Hilfe nach einer beliebigen, der Verteilung zugehörigen Funktion zu suchen oder den Funktionsnamen ohne Klammer auszuführen, z.B.

## function (x, mean = 0, sd = 1, log = FALSE) 
## .Call(C_dnorm, x, mean, sd, log)
## <bytecode: 0x7fb2fe4b2ee8>
## <environment: namespace:stats>

Im Beispiel sehen wir, dass für die Normalverteilung die Parameter mean und sd festgelegt werden müssen und auf die Standardwerte 0 und 1 gesetzt sind. Die Standardwerte werden verwendet, wenn vom Nutzer keine Parameterwerte festgelegt werden.

## [1] 0.2419707
## [1] 0.1994711

2.2 Empirische Verteilungen, Kerneldichte - deskriptive Darstellung

Bei empirischen Daten wollen wir oftmals zunächst die empirische Verteilung unabhängig von einer Verteilungsannahme betrachten. Hierzu bietet es sich an die empirische Dichtefunktion über ein Histogramm bzw. eine Darstellung der nicht-parametrischen Kerneldichteschätzer zu verwenden.

2.3 Satz von Bernoulli, zentraler Hauptsatz der Statistik

Ein wichtiger Aspekt bei der statistischen Modellierung ist der Versuch reale Vorgänge in vereinfachter Form durch Simulationen nachzubilden. Die statistische Rechtfertigung hierfür liegt im Hauptsatz der Statistik bzw. dem Theorem von Bernoulli, welches dem Hauptsatz der Statistik gedanklich vorweg geht. Beim Theorem von Bernoulli geht es darum, dass die relative Häufigkeit eines eines Ereignises bei \(n\)-facher Wiederholung gegen die tatsächliche Wahrscheinlichkeit für das Eintreten des Ereignises konverbiert. Mathematisch ausgedrückt: \[\lim_{n \to \infty} P(\sup | \hat{f}_n(A) - P(A) | < c) = 1, \text{ für beliebige } c \in \mathbb{R}\]

Dies bedeutet, dass wir über relative Häufigkeiten Wahrscheinlichkeiten approximieren können und diese Approximation in der Tendenz genauer werden, wenn wir die Möglichkeit haben ein Ereignis sehr oft zu beobachten. In Form von Simulationen können wir Zufallsvorgänge beliebig oft wiederholen. Ein Beispiel:

## [1] 0.19

Analog lässt sich der Hauptsatz der Statistik durch Simulation zeigen:

2.4 Schätzen von Verteilungen

Wir werden beim schätzen on Verteilungen größtenteils die Maximum-Likelihood-Methode verwenden. Diese ist für die meisten Verteilungen bereits in R bzw. in R Paketen implementiert. Wir wollen es zunächst trotzdem “per Hand” versuchen, um zu zeigen, was in den entsprechenden Funktionen passiert. Geben sei die Situation, dass wir einen fiktiven Datensatz mit 1000 Realsierungen einer Zufallszahl haben (hier z.B. Anzahl der positiven Wertentwicklungen in einer Stunde). Wir gehen davon aus, dass die Daten Poisson-verteilt sind und möchten für diese Daten gemäß der Maximum-Likelihood-Methode den Parameter \(\lambda\) schätzen. Zunächst ziehen wir unseren Datensatz indem wir die Daten als Zufallszahlen ziehen.

##  [1] 4 4 1 3 2 2 3 0 2 3

Wir geben hier den Wert für \(\lambda\) selbst beispielhaft vor. Nun wollen wir die Log-Likelihood dieser Stichprobe bestimmen, gegeben wir setzen Werte für \(\lambda\) ein. Wir schreiben hierfür eine eigene Funktion.

## [1] -1746.267

Nun betrachten wir, wie sich die Log-Likelihood für verschiedene Werte von \(\lambda\) verhält.

Man sieht dass die Log-Likelihood-Funktion ihren Maximalwert in der Nähe des Wertes \(\lambda = 2\) besitzt. Für die Poissonverteilung ist dieser Maximalwert auch rechnerisch ermittelbar. Wir werden jedoch auch oft auf Probleme stoßen, für die keine rechnerische Lösung existiert. Hierfür besteht eine Vielzahl an numerischen Verfahren. Viele dieser Verfahren benötigen die erste Ableitung, sprich den Gradienten, der Log-Likelihood-Funktion. In unserem Fall:

Mit der Information des Gradienten werden stationäre Punkte, also Punkte an denen potentiell ein Optimum existiert, gesucht. Z.B. wandert man bei der Gradientenabstiegsmethode mit Hilfe des Vorzeichens und der Höhe des Gradienten entlang der Funktion. Wir veranschaulichen dies hier an der negativen Log-Likelihood-Funktion, damit wir dem Namen des Verfahren gerecht werden. Experimentieren Sie selbst ein wenig mit der Anzahl der Iterationsschritte und der Learning Rate!

In den seltesten Fällen muss man sich jedoch selbst die Funktion zur Optimierung einer Funktion schreiben. Ganz im Gegenteil, es existieren in R viele Funktionen, die zum Teil unterschiedliche Optimierungsmethoden implementieren. Meist handelt es sich standardisiert um Minimierungen, somit muss bei einer gewünschten Maximierung die ursprüngliche Funktion mit einem negativen Vorzeichen versehen werden, z.B.:

## $minimum
## [1] 1.937999
## 
## $objective
## [1] 1687.768

oder

## $par
## [1] 1.938
## 
## $objective
## [1] 1687.768
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 7
## 
## $evaluations
## function gradient 
##        8       10 
## 
## $message
## [1] "relative convergence (4)"

Für die einfacheren Verteilungen existieren in R jedoch auch Pakete, die direkt eine Maximum-Likelihood-Maximierung ausführen.

##      lambda  
##   1.93800000 
##  (0.04402272)

2.5 AIC, HMI, pp-Plots und qq-Plots

Haben wir eine oder mehrere Verteilungen geschätzt, so wollen wir überprüfen, wie gut sehr die geschätzte Verteilung mit den empirischen Daten im Einklang ist. Dies können wir deskriptiv mit Hilfe von qq- bzw. pp-Plots beurteilen und zusätzlich Kennuzahlen, wie das AIC oder den Harmonic Mass Index berechnen. Im folgenden Beispiel modellieren wir Zufallszahlen, die durch eine Mischung von mehreren Verteilungen zustande kommt, so dass keine Verteilung die “wahre” Verteilung ist. Dann schätzen wir drei Verteilungen an die empirischen Daten.

  • Normalverteilung
  • Student t Verteilung
  • Generalized Error Distribution

Danach vergleichen wir die drei Verteilungen, um uns wie im echten Leben für ein Modell zu entscheiden. Zum Schätzen verwenden wir die ‘fun_fun’-Funktion (s. Code), die Log-Likelihood, AIC, HMI und eine pp- bzw. qq-plot ausgibt. Als Verteilungen können gewählt werden:

  • Normalverteilung: ‘norm’
  • Schiefe Normalverteilung: ‘snorm’
  • Student t Verteilung: ‘std’
  • Schiefe Student t Verteilung: ‘sstd’
  • Generalized Error Disribution: ‘ged’
  • Schiefe Generalized Error Distribution: ‘sged’
#Funktion, welche die angegebene Verteilung schätzt, ll, aic, hmi, pp- und qq-plot ausgibt
library(rugarch)

fun_fun <- function(x, distr = 'norm'){

  distr_set <- c('norm', 'snorm', 'std', 'sstd', 'ged', 'sged')
  
  if(!(any(distr == distr_set))){
    stop(paste('Bitte verwenden Sie eine dieser Verteilungen: norm, snorm, std, sstd, ged, sged'))
  }
  
  
  fit <- fitdist(distribution = distr, x)
  
  
  if(distr == 'norm'){
    
    ll <- sum(log(ddist(distribution = distr, x, mu = fit$pars[1], sigma = fit$pars[2])))
    aic <- -2*ll + 2 * length(fit$pars)
    
      # integrate(function(p1, m, s, data){abs(p1 - pdist(distribution = distr, quantile(data, p1), mu = m, sigma = s))},
                         #lower = 0, upper = 1, m = fit$pars[1], s = fit$pars[2], data = x, subdivisions = 1000)$value
    
    p_theor <- pdist(distribution = distr, sort(x), mu = fit$pars[1], sigma = fit$pars[2])
    p_emp <- rank(sort(x))/length(x)
    
    hmi <- (2 / length(x)) * sum(abs(p_emp - p_theor)) 
    
    pp_plot <- plot_ly() %>%
      add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
      add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    q_theor <- qdist(distribution = distr, rank(sort(x))/(length(x) + 1), mu = fit$pars[1], sigma = fit$pars[2])
    q_emp <- quantile(x, rank(sort(x))/(length(x) + 1))
  
    qq_plot <- plot_ly() %>%
      add_markers(x = q_theor, y = q_emp, name = 'qq-plot', marker = list(color = 'skyblue', symbol = 'x')) %>%
      add_lines(x = c(min(x), max(x)), y = c(min(x), max(x)), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    p <- subplot(pp_plot, qq_plot, margin = 0.05, titleX = TRUE, titleY = TRUE, which_layout = 1)
    
    list_out <- list(ll = ll, aic = aic, hmi = hmi, p = p)
    
  }
  
  if(any(distr == c('std', 'ged'))){
    
    ll <- sum(log(ddist(distribution = distr, x, mu = fit$pars[1], sigma = fit$pars[2], shape = fit$pars[3])))
    aic <- -2*ll + 2 * length(fit$pars)
    #integrate(function(p1, m, s, nu, data){abs(p1 - pdist(distribution = distr, quantile(data, p1), mu = m, sigma = s, shape = nu))},
                #         lower = 0, upper = 1, m = fit$pars[1], s = fit$pars[2], nu = fit$pars[3], data = x, subdivisions = 1000)$value
    
    p_theor <- pdist(distribution = distr, sort(x), mu = fit$pars[1], sigma = fit$pars[2], shape = fit$pars[3])
    p_emp <- rank(sort(x))/length(x)
    
    hmi <- (2 / length(x)) * sum(abs(p_emp - p_theor))
       
    pp_plot <- plot_ly() %>%
      add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
      add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    q_theor <- qdist(distribution = distr, rank(sort(x))/(length(x) + 1), mu = fit$pars[1], sigma = fit$pars[2], shape = fit$pars[3])
    q_emp <- quantile(x, rank(sort(x))/(length(x) + 1))
    
    qq_plot <- plot_ly() %>%
      add_markers(x = q_theor, y = q_emp, name = 'qq-plot', marker = list(color = 'skyblue', symbol = 'x')) %>%
      add_lines(x = c(min(x), max(x)), y = c(min(x), max(x)), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    p <- subplot(pp_plot, qq_plot, margin = 0.05, titleX = TRUE, titleY = TRUE, which_layout = 1)
    
    list_out <- list(ll = ll, aic = aic, hmi = hmi, p = p)
    
  }
  
  if(any(distr == c('snorm', 'sstd', 'sged'))){
    
    ll <- sum(log(ddist(distribution = distr, x, mu = fit$pars[1], sigma = fit$pars[2], skew = fit$pars[3], shape = fit$pars[4])))
    aic <- -2*ll + 2 * length(fit$pars)
     #integrate(function(p1, m, s, skew, nu, data){abs(p1 - pdist(distribution = distr, quantile(data, p1), mu = m, sigma = s, skew = skew, shape = nu))},
                #         lower = 0, upper = 1, m = fit$pars[1], s = fit$pars[2], skew = fit$pars[3], nu = fit$pars[4], data = x, subdivisions = 1000)$value
    
    p_theor <- pdist(distribution = distr, sort(x), mu = fit$pars[1], sigma = fit$pars[2], skew = fit$pars[3], shape = fit$pars[4])
    p_emp <- rank(sort(x))/length(x)
    
    hmi <- (2 / length(x)) * sum(abs(p_emp - p_theor))
    
    pp_plot <- plot_ly() %>%
      add_markers(x = p_theor, y = p_emp, name = 'pp-plot', marker = list(color = 'grey', symbol = 'x')) %>%
      add_lines(x = c(0, 1), y = c(0, 1), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    q_theor <- qdist(distribution = distr, rank(sort(x))/(length(x) + 1), mu = fit$pars[1], sigma = fit$pars[2], skew = fit$pars[3], shape = fit$pars[4])
    q_emp <- quantile(x, rank(sort(x))/(length(x) + 1))
    
    qq_plot <- plot_ly() %>%
      add_markers(x = q_theor, y = q_emp, name = 'qq-plot', marker = list(color = 'skyblue', symbol = 'x')) %>%
      add_lines(x = c(min(x), max(x)), y = c(min(x), max(x)), line = list(color = 'black'), showlegend = FALSE) %>%
      layout(
        xaxis = list(title = 'Theoretische Verteilung'),
        yaxis = list(title = 'Empirische Verteilung')
      )
    
    p <- subplot(pp_plot, qq_plot, margin = 0.05, titleX = TRUE, titleY = TRUE, which_layout = 1)
    
    list_out <- list(ll = ll, aic = aic, hmi = hmi, p = p)
    
  }
  
return(list_out)
  
}

Nun mit den simulierten Daten:

Table 2.1: Übersicht des Datenfits
AIC HMI
Normal 12500.71 0.0153429
Student 12447.76 0.0084073
GED 12466.58 0.0080974

Bei beiden Kennzahlen gilt, je kleiner umso besser, somit kommen wir hier zu keiner eindeutigen Entscheidung, jedoch würde man die Verwendung der Normalverteilung ausschließen. Zusätzlich steht uns noch die grafische Analyse zur Verfügung.

Normalverteilung

Student t Verteilung

Generalized Error Verteilung