6 Induktive Statistik

Eine wahrscheinliche Unmöglichkeit ist immer einer wenig überzeugenden Möglichkeit vorzuziehen.

Aristoteles

Induktive Statistik erlaubt es, statistische Schlüsse mittels stochastischer Modelle zu treffen, um mit Hilfe einer Stichprobe allgemeine Schlussfolgerungen über die Grundgesamtheit ziehen zu können.

6.1 Wahrscheinlichkeitsverteilungen

R hat als ursprüngliches Tool für Statistiker selbstverständlich eine Fülle an statistischen Wahrscheinlichkeitserteilungen für verschiedenste Anwendungen integriert. Auf diese lässt sich sehr leicht zugreifen, um Werte aus der Dichte- bzw. Wahrscheinlichkeitsfunktion oder der Verteilungsfunktion zu bestimmen, Quantile zu berechnen und Zufallszahlen zu erzeugen.

6.1.1 Wiederholung: Zufallsvariablen und deren Wahrscheinlichkeitsverteilungen

Zunächst folgt an dieser Stelle eine kurze Wiederholung der Erkenntnisse aus der Stochastik (aus der Statistikvorlesung), welche im folgenden dann in R genutzt werden können.

Zufallsereignis

In der Stochastik beschäftigen wir uns mit zufälligen Ereignissen. Definition 6.1 beschreibt diese Ereignisse, für die wir anschließend Wahrscheinlichkeiten messen wollen.

Definition 6.1 Die Ausgänge eines Zufallsvorgangs heißen Ergebnisse $\omega_i$, $i=1,\dots,n$. Der Ergebnisraum $\Omega=\{\omega_1,\dots,\omega_n\}$ ist die Menge aller möglichen Ausgänge. Jede Teilmenge von $\Omega$ heißt Zufallsereignis. Die einelementigen Teilmengen von $\Omega$, d. h. $\{\omega_1\}, \dots, \{\omega_n\}$ bezeichnet man als Elementarereignisse.

Wahrscheinlichkeitsmaß

Wollen wir die Wahrscheinlichkeiten von Zufallsereignissen messen, benötigen wir hierzu ein Wahrscheinlichkeitsmaß. Definition 6.2 liefert hierfür die Grundlage.

Definition 6.2 Eine Abbildung $P:\{A:A\in\Omega\}\rightarrow[0,1]$, die allen Ereignissen aus $A\subset\Omega$ eines Zufallsexperimentes eine Zahl $P(A)$ zuordnet und die Bedingungen

$0 \leq P(A) \leq 1$ für alle $A\in\Omega$
$P(\Omega)=1$
$P(A\cup B) = P(A) + P(B)$ falls $A$ und $B$ disjunkt sind, also falls $A\cap B = \emptyset$ gilt

erfüllt heißt Wahrscheinlichkeitsmaß. $P(A)$ nennt man dann Wahrscheinlichkeit von $A$. Die Bedingungen 1-3 heißen auch Kolmogoroff’sche Axiome.

Zufallsvariable

Da Elementarereignisse und deren Wahrscheinlichkeiten häufig zu unflexibel und nicht handhabbar sind, nutzt man in der Regel auf diesen Ereignissen basierende Zufallsvariablen (vgl. Definition 6.3), um mit deren Wahrscheinlichkeiten zu rechnen. Die Unterscheidung zwischen diskreten und stetigen Zufallsvariablen ist dabei sehr wichtig.

Definition 6.3 Eine Funktion $X:\Omega \rightarrow \mathbb{R}$, die jedem Elementarereignis $\omega_i$ eine reelle Zahl $X(\omega_i)$ zuordnet heißt Zufallsvariable.

Die Zahl $X(\omega_i)=x$ heißt Realisation der Zufallsvariable $X$.

Eine Zufallsvariable $X$ heißt diskret, wenn sie endlich viele oder abzählbar viele mögliche Realisationen besitzt.

Eine Zufallsvariable $X$ heißt stetig, wenn ihre möglichen Realisationen alle reellen Zahlen (ggf. innerhalb gewisser Grenzen) sind.

Dichte- bzw. Wahrscheinlichkeitsfunktion

Für diskrete Zufallsvariablen nutzt man die Wahrscheinlichkeitsfunktion, um die Wahrscheinlichkeiten der möglichen Realisationen zu messen (vgl. Definition 6.4)

Definition 6.4 Für eine diskrete Wahrscheinlichkeitsverteilung $X$ mit Realisationen $x_1,x_2,x_3,\dots$ heißt die Funktion \[ f^X(x)=f(x)= \begin{cases} P(X=x_i) & x=x_i\in\{x_1,x_2,x_3,\dots\} \\ 0 & \text{sonst} \end{cases} \] Wahrscheinlichkeitsfunktion von $X$.

Für stetige Zufallsvariablen ist die Wahrscheinlichkeit $P(X=x) = 0$, da es überabzählbar viele mögliche Realisationen gibt. Daher ist die Wahrscheinlichkeitsfunktion für stetige Zufallsvariablen nicht wohl definiert. Stattdessen verwendet man die Wahrscheinlichkeitsdichte aus Definition 6.5, welche statt der Wahrscheinlichkeit für einzelne Realisationen die Wahrscheinlichkeit von Realisationen innerhalb eines Intervalls misst.

Definition 6.5 Eine Zufallsvariable $X$ heißt stetig, wenn es eine Funktion $f(x)\geq 0$ gibt, so dass für jedes Intervall $[a,b]$ \[ P(a\leq X \leq b) = \int_a^b f(x) dx \] gilt. Die Funktion $f(x)$ heißt (Wahrscheinlichkeits-) Dichte von $X$.

Verteilungsfunktion

Die Verteilungsfunktion kumuliert die Wahrscheinlichkeiten auf. Im diskreten Fall geschieht dies wegen der “nur” abzählbar vielen möglichen Realisationen über eine Summe (vgl. Definition 6.6).

Definition 6.6 Für eine diskrete Wahrscheinlichkeitsverteilung $X$ mit Realisationen $x_1,x_2,x_3,\dots$ heißt die Funktion \[ F(x) = P(X\leq x)=\sum_{i: x_i\leq x}^{} f(x_i) \] Verteilungsfunktion der Zufallsvariable $X$.

Bei stetigen Zufallsvariablen ist die Verteilungsfunktion dann wie in Definition 6.7 angegeben als Integral über die Dichtefunktion definiert.

Definition 6.7 Für eine stetige Zufallsvariable $X$ ist die Verteilungsfunktion definiert als \[ F(x) = P(\infty \leq X \leq x) = P(X \leq x) = \int_{-\infty}^x f(t) dt. \]

Die empirische Verteilungsfunktion ist das entsprechende empirische Pendant zur theoretischen Verteilungsfunktion.

Quantile

Die Quantilsfunktion ist eine verallgemeinerte inverse Verteilungsfunktion

Definition 6.8 Für eine Zufallsvariable $X$ ist die verallgemeinerte inverse Verteilungsfunktion definiert als \[ F^{-1}(u) = \inf\{x\in\mathbb{R} | F(x) \geq u\}. \]

Sie ist das entsprechende theoretische Pendant zur Berechnung allgemeiner $p$-Quantile.

6.1.2 Wahrscheinlichkeitsverteilungen in R

In R sind eine Reihe Wahrscheinlichkeitsverteilungen hinterlegt. Tabelle 6.1 respektive Tabelle 6.2 listen die wichtigsten diskreten und stetigen Verteilungsfamilien auf und deren entsprechendes Schlüsselwort auf.

Tabelle 6.1: Diskrete Wahrscheinlichkeitsverteilungen in R

Verteilung	Schlüsslwort
Binomialverteilung	`binom`
Multinomialverteilung	`multinom`
Hypergeometrische Verteilung	`hyper`
Poisson Verteilung	`pois`

Tabelle 6.2: Stetige Wahrscheinlichkeitsverteilungen in R

Verteilung	Schlüsslwort
Gleichverteilung	`unif`
Normalverteilung	`norm`
$t$-Verteilung	`t`
$\chi^2$-Verteilung	`chisq`
$F$-Verteilung	`f`
Log-Normalverteilung	`lnorm`
Exponentialverteilung	`exp`
Cauchy Verteilung	`cauchy`
Multivariate Normal- und $t$-Verteilung	`mvtnorm::pmvnorm`

Dichte- bzw. Wahrscheinlichkeitsfunktion

Um nun die Werte der Wahrscheinlichkeits- bzw. Dichtefunktion zu erhalten, stellt man diesem Schlüsselwort jeweils ein d (engl. density) voran. So erhält man beispielsweise die Wahrscheinlichkeit für einen Wert von $x=1$ aus der Binomialverteilung mit Parametern $n=10$ und $p=0.3$ mit Hilfe von

dbinom(x = 1, size = 10, prob = 0.3)

[1] 0.1210608

Die Wahrscheinlichkeitsfunktion können wir auch grafisch darstellen (vgl. Abbildung 6.1):

x <- 0:10
y <- dbinom(x = x, size = 10, prob = 0.3)
df_binom <- data.frame(x = x, y = y)
ggplot(df_binom, aes(x = x, y = y)) +
  geom_linerange(aes(ymin = 0, ymax = y), linewidth = 1, color = "steelblue4") +
  geom_point(size = 2, color = "steelblue4") +
  labs(x = "x", y = "f(x)")

Abbildung 6.1: Darstellung der Wahrscheinlichkeitsfunktion der Verteilung $ extrm{Bin}(n = 10, p = 0.3)$

Den Wert der Dichtefunktion der Standardnormalverteilung $\phi(x)$ für $x = 0.5$ erhält man über

dnorm(x = 0.5, mean = 0, sd = 1)

[1] 0.3520653

und die grafische Darstellung mit dem folgenden Code, dessen Ergebnis in Abbildung 6.2 dargestellt wird.

x <- seq(-5,5,0.001)
y <- dnorm(x = x)
df_norm <- data.frame(x = x, y = y)
ggplot(df_norm, aes(x = x, y = y)) +
  geom_line(color = "steelblue4", linewidth = 1) +
  labs(x = "x", y = "f(x)")

Abbildung 6.2: Darstellung der Dichte der Standardnormalverteilung

Verteilungsfunktion

Die Werte der Verteilungsfunktion erhält man über ein vorangestelltes p (engl. probability).

So erhält man beispielsweise den Wert der Verteilungsfunktion für $x=5$ aus der Binomialverteilung mit Parametern $n=10$ und $p=0.3$ mit Hilfe von

pbinom(q = 1, size = 10, prob = 0.3)

[1] 0.1493083

Die Verteilungsfunktion können wir auch grafisch darstellen (vgl. Abbildung 6.3):

x <- 0:10
y <- pbinom(q = x, size = 10, prob = 0.3)
df_binom_cdf <- data.frame(x = x, y = y)
ggplot(df_binom_cdf, aes(x = x, y = y)) +
  geom_step(linewidth = 1, color = "steelblue4") +
  geom_point(size = 2, color = "steelblue4") +
  labs(x = "x", y = "F(x)")

Abbildung 6.3: Darstellung der Verteilungsfunktion der Verteilung $ extrm{Bin}(n = 10, p = 0.3)$

Den Wert der Verteilungsfunktion der Standardnormalverteilung $\Phi(x)$ für $x = 0.5$ berechnet man mit

pnorm(q = 0.5, mean = 0, sd = 1)

[1] 0.6914625

und die grafische Darstellung mit dem folgenden Code, dessen Ergebnis in Abbildung 6.4 dargestellt wird.

x <- seq(-5,5,0.001)
y <- pnorm(q = x)
df_norm_cdf <- data.frame(x = x, y = y)
ggplot(df_norm_cdf, aes(x = x, y = y)) +
  geom_line(color = "steelblue4", linewidth = 1) +
  labs(x = "x", y = "F(x)")

Abbildung 6.4: Darstellung der Verteilungsfunktion der Standardnormalverteilung

Quantile

Quantile berechnen wir mit einem vorangestellten q, also z.B. die 97.5%-Quantile für die Binomialverteilung

qbinom(p = 0.975, size = 10, prob = 0.3)

[1] 6

oder für die Standardnormalverteilung

qnorm(p = 0.975, mean = 0, sd = 1)

[1] 1.959964

(Pseudo-)Zufallszahlen

Mit R können wir ebenfalls Pseudozufallszahlen aus den beschriebenen Verteilungen erzeugen. Dazu stellen wir dem Verteilungsnamen ein r voran (engl. random).

So können wir beispielsweise mit

runif(n = 20, min = 0, max = 1)

 [1] 0.4489467747 0.1777305915 0.3258706941 0.1881623317 0.6868427284
 [6] 0.6413182463 0.9649674033 0.4551557782 0.1363010486 0.2547947827
[11] 0.2750000141 0.8956779707 0.1210803553 0.0008099382 0.5835161833
[16] 0.1818776780 0.9147394337 0.2449756043 0.2249650268 0.5305442868

zwanzig gleichverteilte Zufallszahlen aus der Gleichverteilung auf dem Intervall zwischen 0 und 1 erzeugen.

Erzeugen man eine große Zahl standardnormalverteilter Zufallsvariablen und visualisiert die Verteilung durch ein Histogramm, können wir erkennen, dass es sich der theoretischen Dichtefunktion annähert:

x <- rnorm(100000)
df_rnorm <- data.frame(x = x)
ggplot(df_rnorm, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)), bins = 50, fill = "grey80", color = "white") +
  stat_function(fun = dnorm, color = "red", linewidth = 1) +
  labs(x = "x", y = "Dichte")

Abbildung 6.5: Histogramm normalverteilter Zufallsvariablen und Dichtefunktion der Standardnormalverteilung

Tipp

Als Pseudozufall wird bezeichnet, was zufällig erscheint, in Wirklichkeit jedoch berechenbar ist. Pseudozufallszahlen werden auch in R berechnet, d. h. es wird eine Zahlenfolge erzeugt, die zwar zufällig aussieht, es aber nicht ist, da sie durch einen deterministischen Algorithmus berechnet wird. Bei jedem Start der Zufallszahlenberechnung mit gleichem Startwert, der sogenannten Saat (englisch seed), wird die gleiche pseudozufällige Zahlenfolge erzeugt.

Dieser Zusammenhang lässt sich nutzen, um Berechnungen mit Pseudozufallszahlen replizierbar zu machen. Möchte man beispielsweise Ergebnisse von Funktionen wie runif, rnorm, aber auch sample wiederholen und die gleichen Ergebnisse erhalten, muss der Startwert der Zufallszahlenfolge festgelegt werden. Dafür gibt es die Funktion set.seed. Diese Funktion wird zu Anfang des Programms eingefügt und verwendet als einziges Argument eine beliebig zu wählende Zahl. Möchte man das Ergebnis später wiederholen, verwendet man einfach den gleichen “seed”.

6.2 Monte-Carlo-Simulation

Monte-Carlo-Simulationen sind numerische Verfahren, bei denen ein stochastisches Modell sehr oft künstlich wiederholt wird, um Wahrscheinlichkeiten, Erwartungswerte oder Risikomaße näherungsweise zu bestimmen. Ausgangspunkt ist also immer ein Zufallsmechanismus, den wir mit Hilfe von Verteilungen in R beschreiben und dann mit Pseudozufallszahlen realisieren.

6.2.1 Motivation: Warum braucht man Monte Carlo?

Viele interessierende Größen lassen sich zwar klar formulieren, aber nur schwer exakt berechnen. Beispiele sind etwa die Wahrscheinlichkeit kritischer Wartezeiten in einem Logistikzentrum, die zu erwartenden Kosten eines vollständigen Panini-Albums oder ganz allgemein das Risiko extremer Ereignisse in komplexen Prozessen.

Monte-Carlo-Simulationen sind besonders dann nützlich, wenn wir den zugrunde liegenden Zufallsmechanismus modellieren können, aber keine einfache geschlossene Formel für das Ergebnis haben. Statt die gesuchte Wahrscheinlichkeit oder den erwarteten Wert exakt herzuleiten, approximieren wir ihn durch viele künstlich erzeugte Wiederholungen des zugrunde liegenden Zufallsprozesses.

6.2.2 Grundidee: Schätzen durch viele Wiederholungen

Die Grundidee einer Monte-Carlo-Simulation besteht immer aus denselben vier Schritten:

Wir legen Modellannahmen fest.
Wir simulieren eine Realisation des Zufallsprozesses.
Wir berechnen für diesen Durchlauf die interessierende Kennzahl.
Wir wiederholen den Vorgang sehr oft und werten die Ergebnisse gemeinsam aus.

Ein einzelner simulierter Ablauf heißt auch Durchlauf oder Simulation. Erst viele Wiederholungen machen aus diesen Einzelwerten eine brauchbare Schätzung.

Interessiert uns beispielsweise die Wahrscheinlichkeit eines Ereignisses $A$, so können wir in jedem Durchlauf notieren, ob $A$ eingetreten ist oder nicht. Der Anteil der Durchläufe mit Ereigniseintritt ist dann eine empirische Wahrscheinlichkeit bzw. ein Monte-Carlo-Schätzer: \[ \hat{P}(A) = \frac{1}{m}\sum_{r=1}^m I_r, \] wobei $I_r$ im $r$-ten Durchlauf den Wert $1$ annimmt, falls $A$ eintritt, und sonst den Wert $0$.

Interessiert uns dagegen der Erwartungswert einer Größe $Z$, so mitteln wir die simulierten Werte $z_1,\dots,z_m$: \[ \hat{E}(Z) = \frac{1}{m}\sum_{r=1}^m z_r. \] Dieser Mittelwert heißt empirischer Erwartungswert oder ebenfalls Monte-Carlo-Schätzer.

6.2.3 Warum funktioniert das?

Die zentrale Intuition hinter Monte Carlo ist, dass sich bei vielen unabhängigen Wiederholungen relative Häufigkeiten und Mittelwerte stabilisieren. In einfacher Sprache besagt das Gesetz der großen Zahlen, dass der simulierte Durchschnitt bei wachsender Zahl von Wiederholungen typischerweise immer näher an den theoretischen Wert heranrückt.

Monte Carlo liefert also keine exakte Lösung, sondern eine numerische Approximation mit Zufallsfehler. Führt man mehr Wiederholungen durch, werden die Ergebnisse in der Regel stabiler, zugleich steigt aber auch der Rechenaufwand.

Eine sehr kleine Demonstration dieser Idee erhalten wir, wenn wir die Wahrscheinlichkeit für eine Sechs beim fairen Würfel unterschiedlich oft schätzen:

estimate_six <- function(n_rep) {
  mean(sample(1:6, size = n_rep, replace = TRUE) == 6)
}

set.seed(123)
sapply(c(100, 1000, 10000), estimate_six)

[1] 0.2100 0.1660 0.1657

Der theoretische Wert beträgt hier $1/6 \approx 0.167$. Typischerweise sehen wir, dass die Schätzung bei 10.000 Wiederholungen weniger stark schwankt als bei nur 100 Wiederholungen.

6.2.4 Beispiel 1: Panini-Album vervollständigen

Ein sehr anschauliches Beispiel für Monte Carlo ist die Frage, wie teuer es im Mittel ist, ein Panini-Album zu vervollständigen. Für die FIFA World Cup Qatar 2022 Stickerkollektion betrachten wir vereinfachend folgendes Modell:

Das Album enthält 670 verschiedene Sticker.
Eine Tüte enthält 5 Sticker.
Eine Tüte kostet 1,00 Euro.
Es gibt keine Tauschbörse und keine Nachbestellung einzelner Sticker.
Jeder einzelne Sticker einer Tüte wird unabhängig zufällig gezogen; damit sind auch doppelte Sticker innerhalb derselben Tüte möglich.

Ein einzelner Durchlauf besteht nun darin, so lange Tüten zu kaufen, bis alle Sticker mindestens einmal gesammelt wurden. Die Funktion arbeitet dabei in drei Schritten:

Wir starten mit einer leeren Sammlung.
Wir ziehen wiederholt neue Tüten und ergänzen die bisher gesammelten Sticker.
Sobald alle 670 Sticker mindestens einmal vorliegen, brechen wir ab und berechnen die Kosten.

simulate_panini_album <- function(n_stickers = 670,
                                  stickers_per_pack = 5,
                                  price_per_pack = 1) {
  collected <- integer(0)
  n_packs <- 0

  while (length(collected) < n_stickers) {
    # Eine neue Tüte ziehen: Jeder Sticker wird unabhängig gezogen.
    pack <- sample(
      seq_len(n_stickers),
      size = stickers_per_pack,
      replace = TRUE
    )

    # Neue Sticker zur bisherigen Sammlung ergänzen.
    collected <- unique(c(collected, pack))
    n_packs <- n_packs + 1
  }

  # Nach vollständiger Sammlung die wichtigsten Kennzahlen zurückgeben.
  c(
    n_packs = n_packs,
    total_stickers = n_packs * stickers_per_pack,
    total_cost = n_packs * price_per_pack
  )
}

set.seed(2025)
simulate_panini_album()

       n_packs total_stickers     total_cost 
          1019           5095           1019

Ein einzelner Durchlauf ist noch wenig aussagekräftig. Erst viele Wiederholungen geben uns eine Schätzung der zu erwartenden Kosten. Dazu wiederholen wir dieselbe Funktion 5.000-mal und fassen anschließend die simulierten Gesamtkosten zusammen.

set.seed(2025)
# 5.000 unabhängige Album-Simulationen durchführen.
panini_sim <- replicate(5000, simulate_panini_album())
panini_df <- as.data.frame(t(panini_sim))

c(
  Mittelwert = mean(panini_df$total_cost),
  Median = median(panini_df$total_cost),
  `1. Quartil` = quantile(panini_df$total_cost, probs = 0.25),
  `3. Quartil` = quantile(panini_df$total_cost, probs = 0.75)
)

    Mittelwert         Median 1. Quartil.25% 3. Quartil.75% 
      953.5968       924.0000       830.0000      1041.0000

ggplot(panini_df, aes(x = total_cost)) +
  geom_histogram(bins = 40, fill = "grey80", color = "white") +
  labs(x = "Gesamtkosten in Euro", y = "Häufigkeit")

Abbildung 6.6: Verteilung der simulierten Gesamtkosten zur Vervollständigung eines Panini-Albums

Das Ergebnis ist kein exakter Preis, sondern eine empirische Schätzung der zu erwartenden Kosten unter den getroffenen Modellannahmen. Die Verteilung zeigt außerdem, dass die tatsächlichen Kosten von Durchlauf zu Durchlauf spürbar schwanken können.

6.2.5 Beispiel 2: Logistikzentrum mit Wartezeiten

Monte-Carlo-Simulationen lassen sich ebenso nutzen, um realistischere Prozesse mit Ressourcenbeschränkungen zu untersuchen. Betrachten wir dazu ein Logistikzentrum, in dem Transporte zwischen 6:00 Uhr und 18:00 Uhr eintreffen. Die Zwischenankunftszeiten seien exponentialverteilt mit einem Mittelwert von 30 Minuten, die Bearbeitungszeiten ebenfalls exponentialverteilt mit einem Mittelwert von 40 Minuten. Es stehen zwei Mitarbeitende zur Verfügung, und ein Transport gilt als kritisch verspätet, wenn seine Wartezeit mehr als 60 Minuten beträgt.

Ein simulierter Tag kann dann wie folgt erzeugt werden. Auch hier hilft es, die Funktion in Phasen zu lesen:

Zunächst werden alle Ankunftszeiten bis 18:00 Uhr erzeugt.
Danach werden für alle Transporte Bearbeitungszeiten gezogen.
Anschließend wird jeder Transport dem gerade zuerst verfügbaren Mitarbeitenden zugewiesen.
Zum Schluss zählen wir, wie viele Transporte kritisch verspätet waren.

simulate_day <- function(start_hour = 6,
                         end_hour = 18,
                         mean_interarrival = 30,
                         mean_service = 40,
                         n_staff = 2,
                         critical_wait = 60) {
  arrival_times <- numeric(0)
  current_time <- start_hour

  # Ankunftszeiten erzeugen, bis der Planungshorizont endet.
  while (TRUE) {
    current_time <- current_time + rexp(1, rate = 1 / mean_interarrival) / 60

    if (current_time > end_hour) {
      break
    }

    arrival_times <- c(arrival_times, current_time)
  }

  n_transports <- length(arrival_times)

  if (n_transports == 0) {
    return(c(
      n_transports = 0,
      critical_transports = 0,
      critical_day = 0
    ))
  }

  # Für jeden Transport eine Bearbeitungszeit ziehen.
  service_times <- rexp(n_transports, rate = 1 / mean_service) / 60
  available_at <- rep(start_hour, n_staff)
  wait_times <- numeric(n_transports)

  for (i in seq_len(n_transports)) {
    # Immer den Mitarbeitenden wählen, der zuerst wieder frei ist.
    staff_id <- which.min(available_at)
    start_service <- max(arrival_times[i], available_at[staff_id])
    wait_times[i] <- start_service - arrival_times[i]
    available_at[staff_id] <- start_service + service_times[i]
  }

  # Kritische Fälle zählen und einen Tagesindikator zurückgeben.
  critical_transports <- sum(wait_times * 60 > critical_wait)

  c(
    n_transports = n_transports,
    critical_transports = critical_transports,
    critical_day = as.numeric(critical_transports > 0)
  )
}

set.seed(2026)
simulate_day()

       n_transports critical_transports        critical_day 
                 15                   0                   0

Wiederholen wir diese Tagessimulation sehr oft, so erhalten wir einen Monte-Carlo-Schätzer für den Anteil kritischer Tage. Dazu speichern wir pro simuliertem Tag unter anderem, ob mindestens ein kritischer Transport aufgetreten ist.

set.seed(2026)
# 5.000 simulierte Tage erzeugen und in einem Data Frame sammeln.
logistics_sim <- replicate(5000, simulate_day())
logistics_df <- as.data.frame(t(logistics_sim))

c(
  `Mittlere Zahl Transporte` = mean(logistics_df$n_transports),
  `Mittlere Zahl kritischer Transporte` = mean(logistics_df$critical_transports),
  `Anteil kritischer Tage` = mean(logistics_df$critical_day)
)

           Mittlere Zahl Transporte Mittlere Zahl kritischer Transporte 
                            24.0162                              3.3614 
             Anteil kritischer Tage 
                             0.5536

ggplot(logistics_df, aes(x = critical_transports)) +
  geom_histogram(binwidth = 1, boundary = -0.5, fill = "grey80", color = "white") +
  labs(x = "Anzahl kritischer Transporte", y = "Häufigkeit")

Abbildung 6.7: Verteilung der Anzahl kritisch verspäteter Transporte pro simuliertem Tag

Dieses Beispiel zeigt, dass Monte Carlo nicht nur für Sammelprobleme, sondern auch für Prozesse mit Abhängigkeiten und begrenzten Ressourcen eingesetzt werden kann. Gerade wenn analytische Formeln unhandlich werden, ist die Simulation oft ein sehr praktikabler Zugang.

6.2.6 Grenzen und gute Praxis

Monte-Carlo-Ergebnisse hängen immer von den getroffenen Modellannahmen ab.
Schlechte Modellannahmen führen auch bei vielen Wiederholungen zu schlechten Ergebnissen.
Ein einzelner simulierter Durchlauf ist in der Regel nicht aussagekräftig.
Mit set.seed() lassen sich Ergebnisse reproduzierbar machen.
Es ist meist sinnvoll, die Simulationslogik früh in kleine Funktionen zu kapseln.
Berichtet werden sollten Kennzahlen aus vielen Wiederholungen und nicht nur Einzelwerte.

Übung

Verändern Sie im Panini-Beispiel den Tütenpreis oder berücksichtigen Sie zusätzlich ein Starter-Bundle. Wie ändern sich die geschätzten Gesamtkosten?
Vergleichen Sie dieselbe Monte-Carlo-Schätzung für den Würfel mit 100, 1.000 und 10.000 Wiederholungen. Wie verändert sich die Stabilität des Ergebnisses?
Verändern Sie im Logistik-Beispiel die Zahl der Mitarbeitenden oder die mittlere Bearbeitungszeit. Wie wirkt sich dies auf den Anteil kritischer Tage aus?

6.3 Schätzen

6.3.1 Parameterschätzung

Parameterschätzung in R ist denkbar einfach. Kennt man die Schätzfunktion, kann diese in R direkt implementiert werden oder auch die bekannten und existierenden Schätzer verwendet werden wie arithmetisches Mittel (mean), Varianz (var) oder Korrelation (cor)

6.4 Testen

R hat zahlreiche statistische Hypothesentests implementiert. Einige wichtige Standardfälle lassen sich mit wenigen Funktionen direkt abdecken.

6.4.1 Parametertests für Lage- und Streuungsparameter

Die Funktion t.test deckt mehrere wichtige Fälle ab. Für einen Einstichproben-$t$-Test geben wir einen Vektor und den zu testenden Mittelwert mu an:

x <- c(67, 64, 71, 68, 66, 70, 65, 69)
t.test(x, mu = 65)


    One Sample t-test

data:  x
t = 2.8868, df = 7, p-value = 0.02342
alternative hypothesis: true mean is not equal to 65
95 percent confidence interval:
 65.45218 69.54782
sample estimates:
mean of x 
     67.5

Ein einfaches Beispiel für den Mittelwertvergleich von zwei unabhängigen Stichproben ist

t.test(1:10, y = c(7:20))


    Welch Two Sample t-test

data:  1:10 and c(7:20)
t = -5.4349, df = 21.982, p-value = 1.855e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.052802  -4.947198
sample estimates:
mean of x mean of y 
      5.5      13.5

Die Ausgabe auf der Konsole gibt unter anderem den $t$-Wert, die Anzahl Freiheitsgrade, den p-Wert und ein Konfidenzintervall an. Verschiedene Argumente verändern die Durchführung des Tests, wie zum Beispiel die Formulierung der Hypothesen (alternative).

Bei gepaarten Beobachtungen, etwa bei Vorher-Nachher-Messungen, verwenden wir zusätzlich paired = TRUE:

vorher <- c(8.2, 7.9, 6.8, 7.1, 7.5, 6.9, 7.4, 8.1)
nachher <- c(7.8, 7.2, 6.5, 6.8, 7.0, 6.7, 7.1, 7.6)
t.test(vorher, nachher, paired = TRUE)


    Paired t-test

data:  vorher and nachher
t = 7.0553, df = 7, p-value = 0.0002014
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.2659384 0.5340616
sample estimates:
mean difference 
            0.4

Der $F$-Test zum Vergleich zweier Varianzen ist mit der Funktion var.test implementiert. Ein einfaches Beispiel ist

x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var.test(x, y)


    F test to compare two variances

data:  x and y
F = 6.2441, num df = 49, denom df = 29, p-value = 1.345e-06
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
  3.137164 11.747691
sample estimates:
ratio of variances 
          6.244067

6.4.2 Tests für Anteile und kategoriale Daten

Für einen einzelnen Anteil eignet sich der Binomialtest binom.test. Damit lässt sich zum Beispiel prüfen, ob eine Münze mit Wahrscheinlichkeit $p = 0{,}5$ Kopf zeigt:

binom.test(x = 58, n = 100, p = 0.5)


    Exact binomial test

data:  58 and 100
number of successes = 58, number of trials = 100, p-value = 0.1332
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4771192 0.6780145
sample estimates:
probability of success 
                  0.58

Für den Vergleich zweier Erfolgswahrscheinlichkeiten kann prop.test verwendet werden. Das folgende Beispiel vergleicht die Überlebenswahrscheinlichkeiten von Frauen und Männern auf der Titanic:

tab_titanic <- margin.table(Titanic, c("Sex", "Survived"))
erfolge <- tab_titanic[, "Yes"]
n <- rowSums(tab_titanic)

prop.test(x = erfolge, n = n, correct = FALSE)


    2-sample test for equality of proportions without continuity correction

data:  erfolge out of n
X-squared = 456.87, df = 1, p-value < 2.2e-16
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.5643339 -0.4754635
sample estimates:
   prop 1    prop 2 
0.2120162 0.7319149

Für Kontingenztafeln verwenden wir häufig den $\chi^2$-Unabhängigkeitstest mit chisq.test:

x <- matrix(c(12, 5, 7, 7), ncol = 2)
chisq.test(x)


    Pearson's Chi-squared test with Yates' continuity correction

data:  x
X-squared = 0.64112, df = 1, p-value = 0.4233

Die gleiche Funktion kann auch als Anpassungstest verwendet werden, etwa um zu prüfen, ob ein Würfel fair ist:

haeufigkeiten <- c(20, 14, 19, 18, 17, 12)
chisq.test(haeufigkeiten, p = rep(1/6, 6))


    Chi-squared test for given probabilities

data:  haeufigkeiten
X-squared = 2.84, df = 5, p-value = 0.7246

Bei kleinen erwarteten Häufigkeiten in einer $2 \times 2$-Tafel ist fisher.test oft besser geeignet:

tab_small <- matrix(c(1, 9, 11, 3), nrow = 2, byrow = TRUE,
                    dimnames = list(Gruppe = c("Therapie", "Kontrolle"),
                                    Erfolg = c("ja", "nein")))
fisher.test(tab_small)


    Fisher's Exact Test for Count Data

data:  tab_small
p-value = 0.002759
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.0006438284 0.4258840381
sample estimates:
odds ratio 
0.03723312

6.4.3 Tests auf Zusammenhang

Für die Korrelationen aus der Funktion cor können ebenfalls Signifikanztests mit cor.test berechnet werden. Standardmäßig wird die Pearson-Korrelation getestet:

cor.test(mtcars$wt, mtcars$mpg, method = "pearson")


    Pearson's product-moment correlation

data:  mtcars$wt and mtcars$mpg
t = -9.559, df = 30, p-value = 1.294e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9338264 -0.7440872
sample estimates:
       cor 
-0.8676594

Wenn eher ein monotoner als ein linearer Zusammenhang untersucht werden soll, ist ein rangbasierter Test mit method = "spearman" oft sinnvoll:

cor.test(mtcars$wt, mtcars$mpg, method = "spearman", exact = FALSE)


    Spearman's rank correlation rho

data:  mtcars$wt and mtcars$mpg
S = 10292, p-value = 1.488e-11
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
-0.886422

6.4.4 Nichtparametrische Tests

Auch nichtparametrische Tests wie der Wilcoxon-Test (wilcox.test) oder der Kruskal-Wallis-Test (kruskal.test) sind in R implementiert. Für zwei unabhängige Gruppen kann der Wilcoxon-Test beispielsweise so berechnet werden:

tooth_1 <- subset(ToothGrowth, dose == 1)
wilcox.test(len ~ supp, data = tooth_1)

Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
compute exact p-value with ties


    Wilcoxon rank sum test with continuity correction

data:  len by supp
W = 88.5, p-value = 0.00403
alternative hypothesis: true location shift is not equal to 0

Für mehr als zwei Gruppen eignet sich der Kruskal-Wallis-Test:

kruskal.test(len ~ factor(dose), data = ToothGrowth)


    Kruskal-Wallis rank sum test

data:  len by factor(dose)
Kruskal-Wallis chi-squared = 40.669, df = 2, p-value = 1.475e-09