Wahrscheinlichkeitsverteilungen in R

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Wahrscheinlichkeit und Statistik
      • Wahrscheinlichkeitsverteilungen in R

Zum besseren Verständnis werden Kenntnisse über Zufallsvariablen vorausgesetzt, wie sie etwa in Grundbegriffe der Wahrscheinlichkeitsrechnung: Diskrete und stetige Zufallsvariablen vermittelt werden.

Einführung

In den Basis-Paketen von R sind für zahlreiche Wahrscheinlichkeitsverteilungen Funktionen implementiert, die den einfachen Umgang mit den Verteilungen ermöglichen. Diese Funktionen werden hier als Service-Funktionen bezeichnet und es gibt zu jeder Wahrscheinlichkeitsverteilung vier Funktionen, die einer einheitlichen Namenskonvention gehorchen und ähnliche Eingabewerte besitzen. Die Aufgaben der vier Service-Funktionen sind:

1. Berechnung der Wahrscheinlichkeitsdichte f(x).
2. Berechnung der Verteilungsfunktion P(X ≤ x).
3. Berechnung des p-Quantils (Umkehrfunktion der Verteilungsfunktion).
4. Erzeugen einer Stichprobe.

In den folgenden Abschnitten werden vier Wahrscheinlichkeitsverteilungen exemplarisch vorgestellt und dabei typische Anwendungen der Service-Funktionen besprochen. Es handelt sich um zwei diskrete Verteilungen:

• Binomialverteilung
• Poisson-Verteilung

und zwei kontinuierliche Verteilungen:

• kontinuierliche Gleichverteilung
• Normalverteilung.

Nach dem Durcharbeiten der Anwendungs-Beispiele sollte der Leser in der Lage sein, ähnliche Aufgaben mit beliebigen Verteilungen zu bearbeiten.

Die vier Service-Funktionen für Wahrscheinlichkeits-Verteilungen

Im Basis-Paket stats findet man eine große Anzahl von Wahrscheinlichkeits-Verteilungen; nachlesen kann man die komplette Liste unter:

?Distributions

Für jede der Verteilungen werden 4 Funktionen angeboten, die einer einheitlichen Namenskonvention gehorchen und mit denen bestimmte Operationen mit den Verteilungen ausgeführt werden können. Die Erklärungen zu den 4 Funktionen findet man dann nicht unter den 4 Namen der Funktionen, sondern unter dem Namen der Verteilung, also zum Beispiel unter Uniform für die Gleichverteilung oder Poisson für die Poisson-Verteilung. Alle diese Verteilungen sind wie Distributions im Paket stats enthalten.

Für das Beispiel der kontinuierlichen Gleichverteilung (Uniform) sind diese 4 Funktionen und ihre Aufgaben in der folgenden Tabelle beschrieben; dabei steht unif als Abkürzung für uniform, der Anfangsbuchstabe des Funktions-Namens wird in der zweiten Spalte erklärt. Die Eingabewerte (und deren default-Werte) der 4 Funktionen kann man am Listing nach der Tabelle ablesen.

Tabelle 1: Für jede der in Distributions genannten Verteilungen werden 4 Service-Funktionen angeboten, die man an ihrem Anfangsbuchstaben leicht erkennen kann. Die Tabelle zeigt diese 4 Funktionen für die Gleichverteilung auf einem Intervall [min; max], das per default gleich [0; 1] gesetzt ist.

dunif(x, min = 0, max = 1, log = FALSE)
punif(q, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
qunif(p, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
runif(n, min = 0, max = 1)

Mit den Eingabewerten min = 0 und max = 1 wird festgelegt, um welche Gleichverteilung es sich handeln soll. Es wird also ein Intervall für die Wertemenge der Zufallsvariable festgelegt, dessen Werte mit jeweils gleicher Wahrscheinlichkeit angenommen werden; Werte außerhalb des Intervalls werden nicht angenommen (oder mit Wahrscheinlichkeit 0). Der default-Wert für das Intervall ist [0; 1].

Mit den Eingabewerten x und q werden jeweils x-Werte für die Wahrscheinlichkeitsdichte beziehungsweise die Verteilungsfunktion vorgegeben. Beim p-Quantil werden Wahrscheinlichkeiten vorgegeben. Um eine gleich verteilte Stichprobe zu erzeugen, muss man natürlich die Länge n der Stichprobe vorgeben.

Die anderen Eingabewerte sind für erste Anwendungen unerheblich und werden unten nur kurz erklärt.

Die Eingabewerte min = 0 und max = 1 sind natürlich charakteristisch für die Gleichverteilung, bei anderen Verteilungen stehen hier die entsprechenden Eingabewerte, die die Verteilung kennzeichnen.

Die Rückgabewerte der Funktionen werden in folgender Tabelle kurz beschrieben.

Tabelle 2: Rückgabewerte der 4 Funktionen für Wahrscheinlichkeits-Verteilungen.

Aufgrund dieser immer gleichen Struktur, ist es wenig hilfreich, sämtliche Verteilungen ausführlich zu besprechen; hier werden einige ausgewählt:

Tabelle 3: Verteilungen, die in diesem Kapitel besprochen werden, sowie ihr Name, unter dem sie in der R-Dokumentation zu finden sind.

Diskrete Verteilungen

Die Binomialverteilung

Eigenschaften der Binomialverteilung

In Abbildung 1 sind die wichtigsten Eigenschaften der Binomialverteilung zusammengefasst.

Abbildung 1: Die Eigenschaften der Binomialverteilung.

Die typische Anwendung der Binomialverteilung kann man sich an einem Glücksspiel verdeutlichen, das N mal unter gleichen Bedingungen unabhängig voneinander wiederholt wird. Beträgt die Gewinn-Wahrscheinlichkeit bei einem Spiel p, so berechnet sich die Wahrscheinlichkeit dafür, dass man bei N Spielen k mal gewinnt, durch Gleichung (1) in Abbildung 1.

Die Wahrscheinlichkeit bei N Spielen k mal zu gewinnen, wenn die Gewinn-Wahrscheinlichkeit für ein Spiel p beträgt, wird üblicherweise bezeichnet mit:

B(N, p, k).

Aus dem Binomialsatz (Gleichung (2) in Abbildung 1) kann man dann die Normierung der Binomialverteilung (Gleichung (3)) und den Erwartungswert für die Anzahl der gewonnenen Spiele (Gleichung (2) herleiten.

Die vier Service-Funktionen für die Binomialverteilung

Im Paket stats findet man unter Binomial die 4 Service-Funktionen zum Umgang mit der Binomialverteilung:

dbinom(x, size, prob, log = FALSE)
pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)
rbinom(n, size, prob)

• Die Konfiguration der Binomialverteilung geschieht mit den Eingabewerten size und prob. Sie geben die Anzahl der Versuche (oder Spiele) beziehungsweise die Treffer-Wahrscheinlichkeit an. Diese Größen wurden oben immer mit N und p bezeichnet.
• Der Eingabewert x steht für die Anzahl der Treffer (die oben immer mit k bezeichnet wurden). Selbstverständlich kann x ein Vektor sein und dann wird durch dbinom() ein Vektor von Wahrscheinlichkeiten berechnet.
• Mit q wird ebenfalls eine Treffer-Anzahl bezeichnet und q kann wiederum ein Vektor sein. Mit pbinom(q, size, prob) wird dann die Wahrscheinlichkeit berechnet, dass bei size Spielen höchstens q Treffer eintreten.
• Zur Berechnung eines p-Quantils wird mit p eine Wahrscheinlichkeit vorgegeben. Auch p kann ein Vektor sein.
• Das Argument log = FALSE beziehungsweise log.p = FALSE ist per default gleich FALSE gesetzt. Gibt man hier TRUE ein, wird anstelle einer Wahrscheinlichkeit p der natürliche Logarithmus der Wahrscheinlichkeit log(p) berechnet. (Da sich die Wahrscheinlichkeiten oft um mehrere Größenordnungen unterscheiden – siehe etwa Abbildung 2 –, kann dies hilfreich sein, um die Zahlenwerte besser zu vergleichen.)
• Das Argument lower.tail = TRUE gibt an, dass Wahrscheinlichkeiten für die Verteilungsfunktion und p-Quantile als P(X ≤ q) berechnet werden. Setzt man lower.tail = FALSE wird die Wahrscheinlichkeit des Gegenereignisses berechnet.

Im folgenden Unterabschnitt werden einfache Beispiele zur Verwendung der Eingabewerte gezeigt.

Die Bedeutung der Eingabewerte

Wenn man zum ersten Mal mit den Service-Funktionen arbeitet, kann es hilfreich sein, zuerst eine kleine Anwendung zu schreiben, die die Wirkung der Eingabewerte klärt – insbesondere bei diskreten Zufallsvariablen X ist es wichtig zu wissen, ob durch die Verteilungsfunktion (hier pbinom()) eine Wahrscheinlichkeit als

P(X < x) oder P(X ≤ x)

berechnet wird. Im folgenden Skript werden dazu einfache Tests ausgeführt:

# unfaire Münze (Trefferwahrscheinlichkeit p = 0.4):
dbinom(x = c(0, 1), size = 1, prob = 0.4)
# 0.6 0.4

# logarithmische Ausgabe:
dbinom(x = c(0, 1), size = 1, prob = 0.4, log = TRUE)
# -0.5108256 -0.9162907

# 0 oder 1 Treffer mit pbinom():
pbinom(q = 1, size = 1, prob = 0.4)
# 1

# Gegenereignis:
pbinom(q = 1, size = 1, prob = 0.4, lower.tail = FALSE)
# 0

# 0 Treffer:
pbinom(q = 0, size = 1, prob = 0.4)
# 0.6

# Gegenereignis:
pbinom(q = 0, size = 1, prob = 0.4, lower.tail = FALSE)
# 0.4

# p-Quantil:
qbinom(p = 0.6, size = 1, prob = 0.4)
# 0

qbinom(p = 0.7, size = 1, prob = 0.4)
# 1

Man erkennt insbesondere, dass die Funktion pbinom() – einfach konfiguriert – Wahrscheinlichkeiten der Art

P(X ≤ q)

berechnet (siehe Zeile 10).

Aufgabe:

Formulieren Sie entsprechend, wie das p-Quantil qbinom() berechnet wird (siehe Zeile 26 und 29).

Testen Sie den Einsatz von lower.tail = FALSE bei qbinom().

Beispiele zur Anwendung der Binomialverteilung

1. Beispiel: Wahrscheinlichkeiten für die Anzahl der 6 beim sechsmaligen Würfeln

Ein Laplace-Würfel wird 6 mal nacheinander geworfen. Die Würfe seien unabhängig voneinander. Wie groß ist die Wahrscheinlichkeit dafür, dass man k = 0, 1, ..., 6 mal eine 6 würfelt?

Überprüfen Sie, ob die Summe der Wahrscheinlichkeiten tatsächlich 1 ergibt!

Wie groß ist der Erwartungswert der Anzahl der 6 bei 6 Würfen?

Lösung: Die Wahrscheinlichkeiten werden mit Hilfe der Funktion dbinom() berechnet:

dbinom(x, size, prob, log = FALSE)

• Für das Argument x werden die Anzahlen der Treffer eingegeben, also (0:6) .
• Die Anzahl der Würfe wird mit dem Argument size eingegeben.
• Die Treffer-Wahrscheinlichkeit ist die Wahrscheinlichkeit für eine 6, hier also prob = 1/6 .

Die gesuchten Wahrscheinlichkeiten berechnen sich gemäß Zeile 2:

ks <- (0:6)
probs <- dbinom(x = ks, size = 6, prob = 1/6)

print( probs, digits = 3)
# 3.35e-01 4.02e-01 2.01e-01 5.36e-02 8.04e-03 6.43e-04 2.14e-05

# Normierung:
sum(probs)
# 1

# Erwartungswert für die Anzahl der Treffer:
sum(probs * ks)
# 1

In Zeile 8 wird die Summe der Wahrscheinlichkeiten berechnet.

Der Erwartungswert für die Anzahl der 6 wird in Zeile 12 berechnet und ergibt natürlich 1.

Mit den folgenden Anweisungen wird der Plot aus Abbildung 2 erzeugt:

plot(x = ks, y = probs, col = "blue", type = "h", lwd = 3, 
     main = "B(N = 6, p = 1/6, k)", ylab = "P(k)", xlab = "k")
grid()

Abbildung 2: Histogramm für die Wahrscheinlichkeiten, bei 6 Würfen k = 0, 1, ..., 6 mal die 6 zu würfeln. Das Stabdiagramm zeigt die Binomialverteilung mit N = 6 und p = 1/6.

2. Beispiel: Bürgermeisterwahl

Bei der Bürgermeisterwahl in einem Ort mit 1000 Wahlberechtigten gebe es zwei Kandidaten A und B. Es werde angenommen, dass alle Wahlberechtigten an der Wahl teilnehmen. Die Wahl gilt als gewonnenen, wenn ein Kandidat 501 oder mehr Stimmen erhält.

1. Angenommen alle Wahlberechtigten treffen ihre Entscheidung rein zufällig und unabhängig voneinander (etwa durch den Wurf einer fairen Münze). Wie groß ist die Wahrscheinlichkeit dafür, dass A beziehungsweise B die Wahl gewinnt?

2. Nach den Voraussetzungen ist klar, dass der Erwartungswert für den Stimmenanteil jedes Kandidaten bei genau 50 Prozent liegt. Große Abweichungen davon sind sehr unwahrscheinlich. In welchem Intervall (um den Erwartungswert von 50 Prozent) liegt der Stimmenanteil für einen Kandidaten mit einer Wahrscheinlichkeit von 90 Prozent? In welchem Intervall (um den Erwartungswert von 50 Prozent) mit einer Wahrscheinlichkeit von 99 Prozent?

3. Der Kandidat A kann 50 Wahlberechtigte davon überzeugen, ihn zu wählen; die anderen Wahlberechtigten treffen ihre Entscheidung wieder zufällig. Wie groß ist jetzt die Wahrscheinlichkeit dafür, dass A beziehungsweise B gewinnt. (Achtung: es ist hier nicht danach gefragt, welchen Stimmenanteil die Kandidaten erhalten.)

4. Führen Sie (unter den zuletzt genannten Voraussetzungen) solange Simulationen der Wahl aus, bis B zum ersten Mal die Wahl gewinnt oder eine Stichwahl erforderlich ist.

Lösung:

1. Die Wahl wird nur dann nicht entschieden, wenn genau 500 Wahlberechtigte für einen Kandidaten stimmen. Die Wahrscheinlichkeit dafür berechnet sich mit Hilfe von dbinom(), was etwa 0.025 ergibt:

dbinom(x = 500, size = 1000, prob = 0.5)
# 0.02522502

Das Gegenereignis dazu ist, dass die Wahl entschieden wird; aus Symmetrie-Gründen sind die Gewinn-Wahrscheinlichkeiten für beide Kandidaten identisch, daher gewinnt jeder mit einer Wahrscheinlichkeit von etwa 0.487.

p.500 <- dbinom(x = 500, size = 1000, prob = 0.5)
p.500
# 0.02522502
p.A <- (1 - p.500) / 2
p.A
# 0.4873875

Abbildung 3: Histogramm für die möglichen Ergebnisse der Bürgermeisterwahl. Das Diagramm ist folgendermaßen zu lesen: Die Binomialverteilung B(N = 1000, p = 0.5, k) gibt die Wahrscheinlichkeit dafür an, dass ein Kandidat k Stimmen erhält. Die möglichen Werte von k = 0, 1, ..., 1000 sind auf der x-Achse angetragen, die zugehörigen Wahrscheinlichkeiten auf der y-Achse. Das Maximum der Wahrscheinlichkeit liegt bei k = 500 und beträgt etwa 0.0252. Die Bedeutung der roten und grünen Linien wird bei der 2. Frage erklärt.

2. Die Fragestellung, in welchem Intervall sich das Ergebnis mit einer Wahrscheinlichkeit von 90 Prozent befinden wird, kann man leicht an Abbildung 3 veranschaulichen: Die Summe der Wahrscheinlichkeiten (Höhe der Balken) muss 0.9 ergeben – dabei werden die Wahrscheinlichkeiten "in der Mitte des Diagramms" zwischen einem k₁ und k₂ addiert. Und da das Diagramm symmetrisch zu k = 500 ist, muss die Summe der Wahrscheinlichkeiten links von k₁ beziehungsweise rechts von k₂ jeweils 0.05 ergeben. Aber dies ist eine Fragestellung für das p-Quantil: Gesucht ist das p-Quantil der Binomialverteilung zur Wahrscheinlichkeit 0.05. (Entsprechend wählt man die Wahrscheinlichkeit 0.005, wenn das Intervall zu 99 Prozent gesucht ist.)

Hier ist also die Wahrscheinlichkeit 0.05 gegeben und die Schranke k₁ gesucht, so dass die Summe der Wahrscheinlichkeiten der Binomialverteilung unterhalb k₁ gerade 0.05 ergibt. In mathematischer Ausdrucksweise: Gesucht ist das p-Quantil der Binomialverteilung zur Wahrscheinlichkeit 0.05 (bei 1000 Versuchen und Treffer-Wahrscheinlichkeit 0.5). Das p-Quantil wird mit qbinom() berechnet:

qbinom(p = 0.05, size = 1000, prob = 0.5)
# 474

qbinom(p = 0.005, size = 1000, prob = 0.5)
# 459

Das heißt das Wahlergebnis für A liegt mit 90 Prozent Wahrscheinlichkeit im Intervall von 0.475 bis 0.525 und mit 99 Prozent Wahrscheinlichkeit im Intervall von 0.460 bis 0.540. Die entsprechenden Intervalle (in absoluten Zahlen) sind in Abbildung 3 rot beziehungsweise grün gekennzeichnet.

Abbildung 3 wurde erzeugt mit:

N <- 1000
ks <- (0:N)

plot(x = ks, y = dbinom(x = ks, size = N, prob = 0.5), col = "blue", type = "h", 
     xlab = "k", ylab = "P(k)",
     main = "B(N = 1000, p = 0.5, k)", frame.plot = TRUE, pch = 1)
grid()
# 90 Prozent innerhalb:
lines(x = c(475, 475), y = c(0, 0.025), col = "red")
lines(x = c(525, 525), y = c(0, 0.025), col = "red")
# 99 Prozent innerhalb:
lines(x = c(460, 460), y = c(0, 0.025), col = "green")
lines(x = c(540, 540), y = c(0, 0.025), col = "green")

3. Wenn 50 Wähler mit Sicherheit den Kandidaten A wählen, dann reicht es, wenn von den verbleibenden 950 Wahlberechtigten 451 oder mehr den Kandidaten A wählen. Die Wahrscheinlichkeit dafür kann mit pbinom() (Verteilungsfunktion der Binomialverteilung) berechnet werden:

p.A <- 1 - pbinom(q = 450, size = 950, prob = 0.5)
p.A
# 0.9440807

Der Aufruf von pbinom() oben berechnet die Wahrscheinlichkeit dafür, dass von 950 Wahlberechtigten 0, 1, ... oder 450 den Kandidaten A wählen. Gesucht ist aber das Gegenereignis und man erhält etwa 94.4 Prozent Wahrscheinlichkeit dafür, dass A die Wahl gewinnt.

4. Die Simulationen werden mit der Funktion rbinom() in einer repeat-Schleife ausgeführt; rbinom() gibt die Anzahl der Treffer zurück, die jeweils ausgegeben werden (mit den cat()-Befehlen).

repeat{
  z <- rbinom(n = 1, size = 950, prob = 0.5)
  if(z > 450){
    cat("A gewinnt: ", z, "\n")
  } else {
    break()
  }
}
if(identical(z, 450L)){
  cat("Stichwahl erforderlich\n")
} else cat("B gewinnt: ", z, "\n")
# A gewinnt:  488 
# A gewinnt:  493 
# A gewinnt:  457 
# A gewinnt:  469 
# A gewinnt:  492 
# A gewinnt:  475 
# B gewinnt:  448

Der Kandidat A gewinnt die Wahl, wenn er von den 950 Zufallsentscheidungen mehr als 450 Stimmen bekommt. Solange dies der Fall ist, wird die Schleife wiederholt. Ab Zeile 12 ist eine mögliche Realisierung der Simulation gezeigt. Ausgegeben wird dabei immer die Anzahl der Personen, die ihre Entscheidung zufällig treffen und A wählen (hier die Variable z).

Die Poisson-Verteilung

Eigenschaften der Poisson-Verteilung

Man sagt eine Zufallsvariable gehorcht der Poisson-Verteilung, wenn sie die Werte k = 0, 1, 2, ... mit den Wahrscheinlichkeiten in Gleichung (1) in Abbildung 4 annehmen kann; das heißt die Poisson-Verteilung ist vollständig charakterisiert durch den Parameter λ.

Abbildung 4: Die Berechnung der Wahrscheinlichkeiten der Poisson-Verteilung sowie ihre wichtigsten Eigenschaften.

Für k = 0 nimmt die Poisson-Verteilung stets einen echt positiven Wert an (siehe Gleichung (3) in Abbildung 4).

Abbildung 5 zeigt die Wahrscheinlichkeiten der Poisson-Verteilung für verschiedene Werte des Parameters λ.

Abbildung 5: Die Poisson-Verteilung jeweils für k = 0, 1, ..., 10 für 4 verschiedene Werte von λ.

Die Poisson-Verteilung besitzt eine verblüffende Eigenschaft, die man weder der Definition (Gleichung (1) in Abbildung 4) noch den Stabdiagrammen in Abbildung 5 sofort ansieht: Ist eine Zufallsvariable X Poisson-verteilt, dann stimmen sowohl der Erwartungswert E(X) als auch die Varianz Var(X) mit λ überein.

E(X) = λ, Var(X) = λ.

Diese Formeln werden weiter unten hergeleitet (siehe Beispiel 3).

Eine andere bemerkenswerte Eigenschaft betrifft den Zusammenhang mit der Binomialverteilung. Gehorcht eine Zufallsvariable X der Binomialverteilung (zu N Versuchen und Treffer-Wahrscheinlichkeit p), so ist der Erwartungswert von X gleich N · p:

E(X) = N · p.

Für große N sind die Wahrscheinlichkeiten der Binomialverteilung nur mit riesigem Aufwand zu berechnen. Ist aber p sehr klein und setzt man λ = N · p, so kann man die Wahrscheinlichkeiten der Binomialverteilung in sehr guter Näherung durch die Poisson-Verteilung mit Parameter λ berechnen. Die Beispiele 1 und 2 unten sollten Sie unter diesem Aspekt lesen – die exakte Begründung dieser Eigenschaft wird hier nicht gegeben.

Die vier Service-Funktionen für die Poisson-Verteilung

Da die Poisson-Verteilung allein durch den Parameter λ charakterisiert wird, ist dies der einzige Eingabewert in den vier Service-Funktionen, der dazu dient die Verteilung zu konfigurieren. Die anderen Eingabewerte haben die Bedeutung, die schon bei der Binomialverteilung erläutert wurde:

dpois(x, lambda, log = FALSE)
ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)
qpois(p, lambda, lower.tail = TRUE, log.p = FALSE)
rpois(n, lambda)

Beispiele zur Anwendung der Poisson-Verteilung

1. Beispiel: Die Poisson-Verteilung als Näherung der Binomialverteilung

Aufgabe:

Im ersten Beispiel zur Binomialverteilung Wahrscheinlichkeiten für die Anzahl der 6 beim sechsmaligen Würfeln wurden die Wahrscheinlichkeiten für k = 0, 1, 2, ..., 6 Treffer berechnet. In diesem Beispiel ist N · p = 1 und daher könnte nach obiger Behauptung die Poisson-Verteilung als Näherung der Binomialverteilung verwendet werden. Die vorsichtige Ausdrucksweise "könnte verwendet werden" wird gewählt, da nicht klar ist, ob hier N groß genug ist.

1. Berechnen Sie die Wahrscheinlichkeiten für k = 0, 1, 2, ..., 6 mit Hilfe der Poisson-Verteilung (mit Parameter λ = 1) und stellen Sie diese in einem Diagramm den entsprechenden Wahrscheinlichkeiten der Binomialverteilung gegenüber.

2. Berechnen Sie die Summe dieser 7 Wahrscheinlichkeiten der Poisson-Verteilung und den Erwartungswert der Treffer-Anzahl (berechnet mit diesen 7 Wahrscheinlichkeiten). Interpretieren Sie das Ergebnis.

Lösung:

1. Das folgende Skript zeigt die Berechnung der Wahrscheinlichkeiten der Binomial- und Poisson-Verteilung. Die Binomialverteilung wird wie oben mit size = 6, prob = 1/6 konfiguriert (Zeile 2), die Poisson-Verteilung mit lambda = 1 (Zeile 3). In beiden Fällen wird der Vektor xs <- (0:6) für die Treffer-Anzahlen eingesetzt.

xs <- (0:6)
y1 <- dbinom(x = xs, size = 6, prob = 1/6)
y2 <- dpois(x = xs, lambda = 1)

print(y1, digits = 3)
# 3.35e-01 4.02e-01 2.01e-01 5.36e-02 8.04e-03 6.43e-04 2.14e-05

print(y2, digits = 3)
# 0.367879 0.367879 0.183940 0.061313 0.015328 0.003066 0.000511

An den Ausgaben kann man die Unterschiede der Wahrscheinlichkeiten nachvollziehen. Besser gelingt dies mit dem Diagramm, das durch folgendes Skript erzeugt wird (Abbildung 6):

plot(x = xs, y = y1, col = "blue", type = "p", 
     xlab = "k", ylab = "Wahrscheinlichkeit für k Treffer",
     main = "Binomial- und Poisson-Verteilung", frame.plot = TRUE, pch = 1)
grid()

legend("top", legend =c ("Binomial","Poisson"),
       ncol = 2, cex = 0.8, bty = "n",
       col = c("blue","red"), lty = 1, lwd = 2)

points(x = xs, y = y2, col = "red", type = "p")

Im plot()-Befehl werden die Wahrscheinlichkeiten der Binomialverteilung gezeichnet (Zeile 1), mit points() werden die Wahrscheinlichkeiten der Poisson-Verteilung hinzugefügt (Zeile 10).

Abbildung 6: Gegenüberstellung der Wahrscheinlichkeiten, die durch die Binomialverteilung berechnet wurden (blau) und den Näherungen, die mit Hilfe der Poisson-Verteilung berechnet wurden (rot).

2. Normierung und Erwartungswert: Im folgenden Skript werden die oben berechneten Wahrscheinlichkeiten aufsummiert (Zeile 1) und der Erwartungswert der Treffer-Anzahl berechnet (Zeile 5). Für die Binomialverteilung sind die Ergebnisse bekannt und sie beschreiben das Zufallsexperiment richtig (Zeile 9 zeigt nochmals den Erwartungswert).

sum(y2)
# 0.9999168

# Erwartungswert:
sum(y2 * xs)
# 0.9994058

# Erwartungswert der Binomialverteilung:
sum(y1 * xs)
# 1

Für die Poisson-Verteilung ergibt die Summe der Wahrscheinlichkeiten einen Wert, der etwas kleiner ist als 1. Dies rührt daher, dass gemäß der Poisson-Verteilung nicht nur die k-Werte 0, 1, ..., 6 angenommen werden können, sondern beliebige natürliche Zahlen. Erst wenn man alle Wahrscheinlichkeiten summiert, ergibt sich 1. Man erkennt hier, dass die Poisson-Verteilung das Zufallsexperiment nicht richtig beschreibt, da jetzt bei 6 Würfen auch 7, 8, ... Treffer möglich wären. Dies widerspricht aber nicht der Aussage, dass man die Poisson-Verteilung als Näherung der Binomialverteilung einsetzen kann.

Entsprechend ergibt dann auch der Erwartungswert, der mit den Wahrscheinlichkeiten y2 berechnet wird, nicht 1 sondern einen etwas kleineren Wert. Rechnet man mit allen Wahrscheinlichkeiten zu k = 0, 1, ..., erhält man als Erwartungswert der Treffer-Anzahl den Wert des Parameters λ = 1.

2. Beispiel: Anzahl der Personen am Schalter

An einem Schalter treffen pro Stunde im Durchschnitt 100 Personen ein. Wie groß ist die Wahrscheinlichkeit dafür, dass pro Minute 4 oder mehr Personen am Schalter eintreffen?

Führen Sie 1000 Simulationen für die Anzahl der Personen pro Minute aus und berechnen Sie die relativen Häufigkeiten der Ereignisse

"es treffen k Personen pro Minute ein", k = 0, 1, ...

Vergleichen Sie die relativen Häufigkeiten mit den theoretischen Wahrscheinlichkeiten.

Lösung:

Wenn pro 100 im Durchschnitt 100 Personen eintreffen, dann treffen pro Minute 100/60 = 5/3 Personen ein. Man kann die Anzahl der Personen, die pro Minute eintreffen, durch eine Poisson-Verteilung mit λ = 5/3 modellieren.

Das folgende Skript zeigt die beiden Möglichkeiten, wie man die Wahrscheinlichkeit für 4 oder mehr Personen pro Minute berechnen kann:

1 - ppois(q = 3, lambda = 5/3)
# 0.08826715

ppois(q = 3, lambda = 5/3, lower.tail = FALSE)
# 0.08826715

Zeile 1: Das Gegenereignis zu "4 oder mehr Personen pro Minute" sind 0, 1, 2 oder 3 Personen. Mit ppois(q = 3, lambda = 5/3) werden diese 4 Wahrscheinlichkeiten aufsummiert. Die Wahrscheinlichkeit des Gegenereignisses wird durch Subtraktion von 1 gebildet.

Zeile 4: Man kann das Gegenereignis auch direkt mit Hilfe von lower.tail = FALSE berechnen.

Um eine Simulation durchzuführen, wird rpois() mit der Anzahl n = 1000 und dem Parameter λ konfiguriert (Zeile 1). Das Ergebnis ist ein Vektor der Länge 1000; jede Komponente gibt für eine Realisierung an, wie viele Personen in einer Minute eintreffen. Zur Auswertung wird dieser Vektor mit table() ausgegeben (Zeile 2); um relative Häufigkeiten zu erhalten, wird durch 1000 dividiert:

x <- rpois(n = 1000, lambda = 5/3)
table(x) / 1000
# 0     1     2     3     4     5     6     7 
# 0.174 0.314 0.268 0.159 0.060 0.019 0.004 0.002 

print(dpois(x = (0:7), lambda = 5/3), digits = 3)
# 0.1889 0.3148 0.2623 0.1457 0.0607 0.0202 0.0056 0.0013

Zeile 6: Die entsprechenden theoretischen Wahrscheinlichkeiten werden zum Vergleich ausgegeben.

3. Beispiel: Erwartungswert und Varianz der Poisson-Verteilung

Aufgabe: Berechnen Sie für eine Poisson-verteilte Zufallsvariable X mit Parameter λ die Erwartungswerte

E(X), E(X²)

und die Varianz

Var(X).

Führen Sie eine Simulation wie im letzten Beispiel mit n = 1000 und λ = 5/3 durch und vergleichen Sie die empirischen mit den theoretischen Werten E(X), E(X²), Var(X).

Lösung:

Die Berechnung der Erwartungswerte und der Varianz sind in Abbildung 7 ausgeführt.

Abbildung 7: Berechnung des Erwartungswertes und der Varianz der Poisson-Verteilung.

Die Simulation wird im folgenden Skript durchgeführt; es enthält auch die Anweisungen, um den Vergleich zwischen den theoretischen Wahrscheinlichkeiten und den relativen Häufigkeiten graphisch darzustellen (siehe Abbildung 8):

# Simulation: 1000x zählen, wie viele Personen pro min ankommen:
sim.1000 <- rpois(n = 1000, lambda = 5/3)

# in Tabelle verwandeln: Wie oft trifft das Ereignis "X = k" ein?
t.abs <- table(sim.1000)
t.abs
# sim.1000
# 0   1   2   3   4   5   6   7 
# 208 320 243 135  73  16   4   1 

# maximale Anzahl der Personen pro min:
M <- max(sim.1000)
M
# 7

# Tabelle der relativen Häufigkeiten:
t.rel <- t.abs / 1000
t.rel
# sim.1000
#     0     1     2     3     4     5     6     7 
# 0.208 0.320 0.243 0.135 0.073 0.016 0.004 0.001 

# Plot: empirische Werte als Stabdiagramm,
# theoretische Wahrscheinlichkeit als Messpunkt:
plot(x = t.rel, col = rainbow(M+1), lwd = 3, xlab = "k", ylab = "P(X = k)", main = "Anzahl der Personen pro Minute")
grid()
legend("topright", legend =c ("empirisch","theoretisch"),
       ncol = 2, cex = 0.8, bty = "n", lty = 1, lwd = 2, pch = c(3,1))
points(x = (0:(M+1)), y = dpois(x = (0:(M+1)), lambda = 5/3), 
       col = rainbow(M+1), type = "p", lwd = 2)

# Auswertung der Stichprobe:
# Erwartungswert
E <- sum((0:M) * t.rel) 
E
# 1.614

# E(X^2)
E2 <-  sum(((0:M)*(0:M)) * t.rel) 
E2
# 4.268

# Varianz:
V <- E2 - E * E
V
# 1.663004

Abbildung 8: Die aus der Stichprobe berechneten relativen Häufigkeiten sind als Stabdiagramm dargestellt, die theoretischen Wahrscheinlichkeiten als Messpunkte. Da in der Stichprobe maximal 7 Personen pro Minute am Schalter eintreffen, endet die Darstellung bei k = 7; die theoretischen Wahrscheinlichkeiten sind auch für k > 7 von 0 verschieden.

Kontinuierliche Verteilungen

Die kontinuierliche Gleichverteilung

Definition der kontinuierlichen Gleichverteilung

Paradebeispiele für diskrete Gleichverteilungen sind

1. Das Werfen einer fairen Münze.
2. Das Würfeln mit einem Laplace-Würfel.

Im ersten Fall erscheinen die beiden Elementar-Ereignisse Kopf beziehungsweise Zahl jeweils mit Wahrscheinlichkeit 1/2; im zweiten Fall besitzt jedes der Elementar-Ereignisse 1, 2, ..., 6 die Wahrscheinlichkeit 1/6.

Die kontinuierliche Gleichverteilung kann man sich etwa beim Schießen auf eine Zielscheibe veranschaulichen: Geht man davon aus, dass

1. die x-Koordinate des Treffers nur in einem Intervall [min; max] liegen kann und
2. jede dieser x-Koordinaten mit gleicher Wahrscheinlichkeit angenommen wird,

so kann man die zugehörige Wahrscheinlichkeitsverteilung nicht mehr mit einer diskreten Gleichverteilung modellieren, da jetzt überabzählbar viele Elementar-Ereignisse vorliegen.

Die Wahrscheinlichkeiten von Ereignissen der Art

a ≤ x ≤ b

werden jetzt mit Hilfe einer Wahrscheinlichkeitsdichte f(x) definiert, indem man die Wahrscheinlichkeitsdichte über das Intervall [a; b] integriert.

Im Fall der Gleichverteilung auf dem Intervall [min; max] lautet die Wahrscheinlichkeitsdichte:

f(x) = 1_{[min; max]} (x) / (max - min).

Die Indikatorfunktion 1_{[min; max]} (x) des Intervalls [min; max] nimmt die Werte 1 oder 0 an, je nachdem ob x im Intervall [min; max] liegt oder nicht. Der Faktor 1 / (max - min) sorgt dafür, dass die Wahrscheinlichkeitsdichte richtig normiert ist: bei Intergration über alle x-Werte ergibt sich 1.

Die kontinuierliche Gleichverteilung in R

Die Beschreibung der vier Service-Funktionen zur kontinuierlichen Gleichverteilung findet man im Paket stats unter Uniform:

dunif(x, min = 0, max = 1, log = FALSE)
punif(q, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
qunif(p, min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
runif(n, min = 0, max = 1)

Die Eingabewerte müssen nicht nochmals erklärt werden: Zur Konfiguration der Gleichverteilung reicht es, die Intervallgrenzen min und max anzugeben, alle anderen Eingabewerte haben die Bedeutung, die bereits erläutert wurde.

Anwendungsbeispiel für die kontinuierliche Gleichverteilung

Aufgabe: Stellen Sie die Wahrscheinlichkeitsdichten für die Gleichverteilungen auf den Intervallen [0; 1] beziehungsweise [0; 2] nebeneinander in einem Diagramm dar.

Ebenso die zugehörigen Verteilungsfunktionen.

In einer Simulation soll für diese Gleichverteilungen jeweils eine Stichprobe der Länge 1000 gezogen werden. Für die Stichproben sollen Mittelwert, Varianz und Standardabweichung berechnet werden.

Lösung:

Das folgende Skript erzeugt die beiden Wahrscheinlichkeitsdichten, die in Abbildung 9 nebeneinander dargestellt werden:

par(mfrow = c(1, 2))
xs <- c(-1, 3)
ys <- c(-0.2, 1.2)

curve(expr = dunif(x), xlim = xs, ylim = ys, n = 1001,
      col = "blue", main = "Gleichverteilung auf [0; 1]")
grid()

curve(expr = dunif(x, max = 2), xlim = xs, ylim = ys, n = 1001,
      col = "red", main = "Gleichverteilung auf [0; 2]")
grid()

Zeile 1: Die beiden Plots sollen nebeneinander dargestellt werden, dazu wird mfrow mit einer Zeile und zwei Spalten konfiguriert.

Zeile 2 und 3: Festlegen der x- und y-Bereiche, in denen die Wahrscheinlichkeitsdichten gezeichnet werden sollen.

Zeile 5 und 6: An die Funktion curve() wird mit expr die Wahrscheinlichkeitsdichte der Gleichverteilung auf [0; 1] übergeben. Dazu muss dunif() nicht weiter konfiguriert werden, da min = 0 und max = 1 die default-Werte sind. Die Anzahl der Stützstellen zum Zeichnen wird mit n = 1001 festgelegt; wählt man n zu klein, ist der Anstieg der Flanke in der Indikatorfunktion nicht steil genug.

Zeile 9 und 10: Für die Gleichverteilung auf [0; 2] muss lediglich max = 2 in dunif() gesetzt werden.

Abbildung 9: Die Wahrscheinlichkeitsdichte der Gleichverteilung für die Intervalle [0; 1] und [0; 2].

Zur Darstellung der Verteilungsfunktionen wird lediglich dunif() durch punif() ersetzt (siehe Abbildung 10):

par(mfrow = c(1, 2))

curve(expr = punif(x), xlim = xs, ylim = ys, n = 1001,
      col = "blue", main = "Verteilungsfunktion zur Gleichverteilung auf [0; 1]")
grid()

curve(expr = punif(x, max = 2), xlim = xs, ylim = ys, n = 1001,
      col = "red", main = "Verteilungsfunktion zur Gleichverteilung auf [0; 2]")
grid()

Abbildung 10: Die Verteilungsfunktionen der Gleichverteilung für die Intervalle [0; 1] beziehungsweise [0; 2].

Im folgende Skript werden Stichproben der Länge 1000 zu den beiden besprochenen Gleichverteilungen erzeugt und ausgewertet:

N <- 1000

# Gleichverteilung auf [0; 1]:
v <- runif(n = N)
mean(v)
# 0.4882555
var(v)
# 0.08493159
sd(v)
# 0.2914302

summary(v)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.0002389 0.2294261 0.4803411 0.4882555 0.7439418 0.9984908 

# Gleichverteilung auf [0; 2]:
v <- runif(n = N, max = 2)
mean(v)
# 0.9940314
var(v)
# 3507778
sd(v)
# 0.592265

summary(v)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.000747 0.477173 1.004906 0.994031 1.512672 1.999526

Die Normalverteilung

Die Definition der Normalverteilung

Mit Normalverteilung wird eine Klasse von Wahrscheinlichkeitsverteilungen bezeichnet, deren Wahrscheinlichkeitsdichten die Form einer Gaußschen Glockenkurve besitzen (siehe Abbildung 13). Sie unterscheiden sich in den beiden Parametern Erwartungswert μ und Standardabweichung σ. Speziell bei der Standard-Normalverteilung ist μ = 0 und σ = 1.

Die wichtigsten Eigenschaften der Standard-Normalverteilung sind in Abbildung 11 zusammengefasst:

1. Die Wahrscheinlichkeitsdichte in Gleichung 1.
2. Die Normierung der Wahrscheinlichkeitsdichte in Gleichung 2.
3. Die formale Definition der Verteilungsfunktion in Gleichung 3 (man beachte, dass man hier keine Stammfunktion angeben kann, die aus elementaren Funktionen aufgebaut ist).
4. Die Symmetrie-Eigenschaften der Wahrscheinlichkeitsdichte und der Verteilungsfunktion in Gleichung 4.

Abbildung 11: Definition und Eigenschaften der Standard-Normalverteilung.

Die Normalverteilung mit den Parametern μ und σ wird meist mit N(μ, σ) bezeichnet. Ihre Wahrscheinlichkeitsdichte ist in Gleichung (1) in Abbildung 12 gezeigt.

Der Zusammenhang zwischen einer beliebigen Normalverteilung N(μ, σ) und der Standard-Normalverteilung N(0, 1) wird durch die sogenannte Standardisierung einer Zufallsvariable hergestellt (siehe Gleichung (2) in Abbildung 12).

Abbildung 12: Eigenschaften der Normalverteilung N(μ, σ) und ihr Zusammenhang mit der Standard-Normalverteilung.

Man beachte, dass durch die Bezeichnung μ und σ der Parameter in der Wahrscheinlichkeitsdichte der Normalverteilung deren inhaltliche Bedeutung als Erwartungswert und Standardabweichung lediglich suggeriert wird. Dass sie tatsächlich diese Bedeutung haben, muss aber erst gezeigt werden.

Aufgabe: Zeigen Sie, dass für eine normal verteilte Zufallsvariable X mit der Wahrscheinlichkeitsdichte aus Gleichung (1) in Abbildung 12 gilt:

E(X) = μ, Var(X) = σ².

Die vier Service-Funktionen für die Normalverteilung

Die Beschreibung der vier Service-Funktionen der Normalverteilung findet man im Paket stats unter Normal. Die Eingabewerte der Service-Funktionen lauten:

dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)

Man erkennt, dass die beiden Eingabewerte mean und sd verwendet werden, um die Verteilung zu charakterisieren (es handelt sich dabei um den Mittelwert μ und die Standardabweichung σ der Verteilung); die anderen Eingabewerte müssen nicht mehr erklärt werden.

Aufgabe: Testen Sie das Verhalten der vier Service-Funktionen, wenn sd = 0 gesetzt wird.

Anwendungs-Beispiele

1. Beispiel:

Erzeugen Sie ein Diagramm, in dem nebeneinander dargestellt sind:

1. die Wahrscheinlichkeitsdichte der Standard-Normalverteilung N(0, 1),
2. die Wahrscheinlichkeitsdichte der Normalverteilung N(2, 0.5), also mit Mittelwert μ = 2 und Standardabweichung σ = 0.5.

Erzeugen Sie ein weiteres Diagramm, das die entsprechenden Verteilungsfunktionen zeigt.

Führen Sie zu den beiden Verteilungen Simulationen durch:

• Es sollen Stichproben der Länge n = 1000 erzeugt werden.
• Die Stichproben sollen ausgewertet werden, indem gezählt wird, wie oft Zahlen in den Intervallen [-4.5; -3.5], ..., [3.5; 4.5] vorkommen.

Erzeugen Sie zu der Auswertung ein geeignetes Diagramm.

Lösung:

Das folgende Skript erzeugt das Diagramm für die beiden Wahrscheinlichkeitsdichten (siehe Abbildung 13).

par(mfrow = c(1, 2))
xs <- c(-3, 5)
ys <- c(-0.1, 0.9)

curve(expr = dnorm(x), xlim = xs, ylim = ys, n = 1001,
      col = "red", main = "Standardnormalverteilung N(0, 1)")
grid()

curve(expr = dnorm(x, mean = 2, sd = 0.5), xlim = xs, ylim = ys, n = 1001,
      col = "blue", main = "Normalverteilung N(2, 0.5)")
grid()

Zeile 5: Für die Standard-Normalverteilung N(0, 1) muss die Funktion dnorm() nicht konfiguriert werden, da Mittelwert μ = 0 und Standardabweichung σ = 1 die default-Werte sind.

Zeile 9: Für die Wahrscheinlichkeitsdichte von N(2, 0.5) müssen die Argumente mean und sd in dnorm() gesetzt werden.

Abbildung 13: Die Wahrscheinlichkeitsdichten der Standard-Normalverteilung N(0, 1) und der Normalverteilung N(2, 0.5).

Die zugehörigen Verteilungsfunktionen werden im folgenden Skript erzeugt; zum besseren Vergleich ist die Gerade y = 0.5 eingezeichnet (Zeile 6 und 11):

xs <- c(-3, 5)
ys <- c(-0.1, 1.1)

curve(expr = pnorm(x), xlim = xs, ylim = ys, n = 1001,
      col = "red", main = "Verteilungsfunktion von N(0, 1)")
abline(a = 0.5, b = 0, col = "gray")
grid()

curve(expr = pnorm(x, mean = 2, sd = 0.5), xlim = xs, ylim = ys, n = 1001,
      col = "blue", main = "Verteilungsfunktion von N(2, 0.5)")
abline(a = 0.5, b = 0, col = "gray")
grid()

Im Vergleich zum Skript für die Wahrscheinlichkeitsdichten muss nur dnorm() durch pnorm() ersetzt werden.

Abbildung 14: Die Verteilungsfunktionen der Standard-Normalverteilung N(0, 1) und der Normalverteilung N(2, 0.5).

Die Stichprobe und ihre Auswertung wird im folgenden Skript erzeugt:

brk <- (-5:4) + 0.5

v.0 <- rnorm(n = 1000)
t.v.0 <- table(cut(x = v.0, breaks = brk))

v.2 <- rnorm(n = 1000, mean = 2, sd = 0.5)
t.v.2 <- table(cut(x = v.2, breaks = brk))

par(mfrow = c(1, 2))
ys <- c(0, 800)

plot(t.v.0, col = "red", main = "Standardnormalverteilung N(0, 1)", 
     lwd = 3, ylim = ys)
plot(t.v.2, col = "blue", main = "Normalverteilung N(2, 0.5)",
     lwd = 3, ylim = ys)

Zeile 1: Mit brk werden die Intervallgrenzen definiert, für die später die Histogramme gezeichnet werden.

Zeile 3 und 6: Die Stichprobe wird jeweils mit rnorm() erzeugt; dabei entsteht ein Vektor der Länge 1000.

Zeile 4 und 7: Der Vektor wird in einen Faktor und sofort in eine Tabelle umgewandelt. Hier wird als Auswertung lediglich ein Histogramm erzeugt, daher muss man den Faktor nicht abspeichern. Für weitere Auswertungen würde man ihn verwenden. (Wie man Faktoren mit Hilfe von cut() erzeugt und dabei mit dem Argument breaks geeignete Faktor-levels erzeugt, wird in Faktoren in R: Anwendungen erklärt.)

Abbildung 15: Darstellung der absoluten Häufigkeiten der Stichprobenwerte in den durch brk definierten Intervallen.

Zum besseren Vergleich werden in der folgenden Abbildung die Stichproben und die Wahrscheinlichkeitsdichten zusammen dargestellt.

Abbildung 16: Die absoluten Häufigkeiten werden in relative Häufigkeiten umgerechnet und zusammen mit den Wahrscheinlichkeitsdichten aufgetragen.

Dass in der rechten Abbildung so deutliche Unterschiede zwischen den theoretischen Wahrscheinlichkeitsdichten und den relativen Häufigkeiten bestehen, liegt an der Größe der Intervalle (Länge gleich 1, siehe Abbildung 14 oder die Definition von brk), mit denen die Faktoren gebildet wurden; um eine bessere Annäherung zu erreichen, müsste man die Intervalle deutlich kleiner machen.

Abbildung 16 wurde mit folgendem Skript erstellt:

# Tabellen t.v.0 und t.v.2 wie oben
par(mfrow = c(1, 2))
xs <- c(-4, 4)
ys <- c(-0.2, 0.8)

curve(expr = dnorm(x), xlim = xs, ylim = ys, n = 1001,
      col = "red", main = "Standardnormalverteilung N(0, 1)")
grid()
points(x = -4:4, y = t.v.0 / 1000, type = "h", lwd = 3, col = "red")

curve(expr = dnorm(x, mean = 2, sd = 0.5), xlim = xs, ylim = ys, n = 1001,
      col = "blue", main = "Normalverteilung N(2, 0.5)")
grid()
points(x = -4:4, y = t.v.2 / 1000, type = "h", lwd = 3, col = "blue")

Mit curve() werden wieder die Wahrscheinlichkeitsdichten gezeichnet (Zeile 6 und 11).

Damit man die Ergebnisse der Stichproben mit den Wahrscheinlichkeitsdichten vergleichen kann, muss man die absoluten Häufigkeiten, die in den Tabellen gepeichert sind, in relative Häufigkeiten umrechnen (Division durch 1000 in Zeile 9 und 14).

2. Beispiel: Berechnung der Wahrscheinlichkeiten innerhalb der Standardabweichung, der doppelten beziehungsweise dreifachen Standardabweichung

Wie groß sind für eine gemäß N(μ, σ) verteilte Zufallsvariable X die Wahrscheinlichkeiten:

P(X ∈ [μ - σ; μ + σ]),

P(X ∈ [μ - 2 · σ; μ + 2 · σ]),

P(X ∈ [μ - 3 · σ; μ + 3 · σ])?

Lösung:

Da alle Normalverteilung durch Standardisierung auf die Standard-Normalverteilung zurückgeführt werden können, reicht es diese Frage für die Standard-Normalverteilung zu beantworten. Die entsprechenden Integrale kann man aber nicht analytisch berechnen, daher bleibt nur die Möglichkeit, die Service-Funktionen der Normalverteilung einzusetzen:

pnorm(q = 1) - pnorm(q = -1)
# 0.6826895
pnorm(q = 2) - pnorm(q = -2)
# 0.9544997
pnorm(q = 3) - pnorm(q = -3)
# 0.9973002

Die gesuchten Wahrscheinlichkeiten betragen etwa 68.2, 95.4 und 99.7 Prozent.

Aufgabe: Veranschaulichen Sie sich diese Ergebnisse an den Abbildungen 13 und 14!