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.64485362.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.19947112.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.
#Lade plotly Paket zu Darstellung des Histogramms
library(plotly)
#Ziehe 1000 standardnormalverteilte Zufallszahlen
set.seed(42)
x <- rnorm(1000)
#Schätzen der Kerneldichten mit Standardnormalverteilungskernel und Bandweite 0.1
kd01 <- density(x, kernel = 'gaussian', bw = 0.1)
#Schätzen der Kerneldichten mit Standardnormalverteilungskernel und Bandweite 0.2
kd02 <- density(x, kernel = 'gaussian', bw = 0.2)
#Schätzen der Kerneldichten mit Standardnormalverteilungskernel und Bandweite 0.5
kd05 <- density(x, kernel = 'gaussian', bw = 0.5)
plot_ly() %>%
add_histogram(x = x,
histnorm = 'probability density',
xbins = list(start = -5, end = 5, size = 0.5),
name = 'Histogramm',
marker = list(color = 'lightgrey')
) %>%
add_lines(x = kd01$x, y = kd01$y, name = 'bw = 0.1', line = list(color = 'olivedrab')) %>%
add_lines(x = kd02$x, y = kd02$y, name = 'bw = 0.2', line = list(color = 'purple')) %>%
add_lines(x = kd05$x, y = kd05$y, name = 'bw = 0.5', line = list(color = 'skyblue')) %>%
layout(
title = 'Empirische Verteilung',
xaxis = list(title = TeX('$x$')),
yaxis = list(title = TeX('$ \\hat{f}(x) $'))
) %>%
config(last_plot(), mathjax = 'local')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:
#Die reale Wahrscheinlichkeit, dass eine Standardnormalverteilte Zufallszahl zwischen 0 und 0.5 liegt ist:
p_real <- pnorm(0.5) - pnorm(0)
round(p_real, 2)## [1] 0.19#Anzahl der Wiederholung eines Zufallsvorgangs, das bedeutet wir ziehen n Zufallszahlen und beobachten, wie viele im Intervall [0, 0.5] liegen:
n <- c(10, 20, 50, 100, 500, 1000, 5000, 10000, 50000, 100000)
freq_p <- c()
set.seed(42)
for(i in 1:length(n)){
x_temp <- rnorm(n[i])
freq_p[i] <- mean(x_temp <= 0.5 & x_temp >= 0)
}
plot_ly() %>%
add_lines(x = n, y = freq_p, line = list(color = 'grey'), name = 'rel. Häufigkeit') %>%
add_markers(x = n, y = freq_p, marker = list(color = 'grey'), showlegend = FALSE) %>%
add_lines(x = c(0, 100000), y = c(p_real, p_real), line = list(color = 'olivedrab', dash = 'dash'), name = 'Wahrscheinlichkeit')Analog lässt sich der Hauptsatz der Statistik durch Simulation zeigen:
#ziehe 1000 Zufallszahlen
set.seed(42)
x <- rnorm(1000)
#Wir betrachten die rel. Häufigkeiten, das x <= x_p
x_p <- seq(-5, 5, 0.5)
x_freq <- c()
for(i in 1:length(x_p)){
x_freq[i] <- mean(x <= x_p[i])
}
plot_ly()%>%
add_lines(x_p, x_freq, line = list(color = 'grey'), name = 'Empirische Verteilungsfunktion') %>%
add_lines(x_p, pnorm(x_p), line = list(color = 'olivedrab', dash = 'dash'), name = 'Kum. Standardormalverteilung')%>%
layout(
xaxis = list(title = TeX('$x$')),
yaxis = list(title = TeX('$\\hat{F}(x),~F(x)$'))
) %>%
config(last_plot(), mathjax = 'local')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 3Wir 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.
ll <- function(data, theta = 1){
#Dies sichert dass eine andere Person nur zulässige Werte für lambda eingeben würde
#Sie können das aber im Moment ignorieren;)
if(theta <= 0){stop('Der Wert für lambda muss eine positive Zahl sein')}
log_l <- sum(log(dpois(data, lambda = theta)))
return(log_l)
}
#z.B. ist der Wert der Log-Likelihood für lambda = 1.5
ll(x, theta = 1.5)## [1] -1746.267Nun betrachten wir, wie sich die Log-Likelihood für verschiedene Werte von \(\lambda\) verhält.
theta <- seq(0.1, 3, 0.05)
log_out <- c()
for(i in 1:length(theta)){
log_out[i] <- ll(x, theta[i])
}
plot_ly() %>%
add_lines(x = theta, y = log_out, line = list(color = 'olivedrab')) %>%
layout(
xaxis = list(title = TeX('$\\lambda$')),
yaxis = list(title = TeX('$\\sum_{i = 1}^n \\ln \\left( f(x_i | \\lambda) \\right) $'))
) %>%
config(last_plot, mathjax = 'local')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!
#Startwert
theta_0 <- 0.5
#zugehöriger Wert der negativen Log-Likelihood
ll_neg_grad <- -ll(x, theta_0)
#Lernrate
eta <- 0.0001
#Anzahl der Iterationsschritte
iter <- 20
for(i in 1:iter){
theta_0 <- c(theta_0, theta_0[i] - eta * (-grad_ll(data = x, theta = theta_0[i])))
ll_neg_grad <- c(ll_neg_grad, -ll(x, theta_0[i + 1]))
}
plot_ly() %>%
add_lines(x = theta, y = -log_out, line = list(color = 'olivedrab'), name = 'negative LL') %>%
add_markers(x = theta_0, y = ll_neg_grad, marker = list(color = 'purple'), name = 'Iterationsschritte') %>%
layout(
xaxis = list(title = TeX('$\\lambda$')),
yaxis = list(title = TeX('$-\\sum_{i = 1}^n \\ln \\left( f(x_i | \\lambda) \\right) $'))
) %>%
config(last_plot, mathjax = 'local')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.:
neg_ll <- function(data, theta){-ll(data, theta)}
optimize(neg_ll, interval = c(exp(-10), 5), data = x)## $minimum
## [1] 1.937999
##
## $objective
## [1] 1687.768oder
## $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:
#Simulation der Daten, wir simulieren jeweils 1000 Datenpunkte aus einer der drei Verteilungen
set.seed(42)
x <- c(rdist('norm', n = 1000, mu = 1, sigma = 2), rdist('std', n = 1000, mu = 1, sigma = 2, shape = 3), c(rdist('ged', n = 1000, mu = 1, sigma = 2, shape = 4)))
#Vorab können wir die Verteilung visualisieren
#Kerneldichteschätzer
kd <- density(x, kernel = 'gaussian')
plot_ly() %>%
add_histogram(x = x,
histnorm = 'probability density',
xbins = list(start = -5, end = 5, size = 0.5),
name = 'Histogramm',
marker = list(color = 'lightgrey')
) %>%
add_lines(x = kd$x, y = kd$y, name = 'Kerneldichteschätzer', line = list(color = 'olivedrab')) %>%
layout(
title = 'Empirische Verteilung',
xaxis = list(title = TeX('$x$')),
yaxis = list(title = TeX('$ \\hat{f}(x) $'))
) %>%
config(last_plot(), mathjax = 'local')#Schätze Normalverteilung
norm_fit <- fun_fun(x, distr = 'norm')
#Schätze Studentverteilung
student_fit <- fun_fun(x, distr = 'std')
#Schätze GED-Verteilung
ged_fit <- fun_fun(x, distr = 'ged')
df <- data.frame('AIC' = c(norm_fit$aic, student_fit$aic, ged_fit$aic),
'HMI' = c(norm_fit$hmi, student_fit$hmi, ged_fit$hmi),
row.names = c('Normal', 'Student', 'GED'))
library(knitr)
kable(df, caption = 'Ü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