Hinweis: Die Gliederung zu den Befehlen Und Ergebnissen in R orientiert sich an dem Inhaltsverzeichnis der 17. Auflage der ‘Angewandten Statistik’. Nähere Hinweise zu den verwendeten Formeln sowie Erklärungen zu den Beispielen sind in dem Buch (Ebook) nachzulesen!

zur Druckversion

Hinweis: Die thematische Gliederung zu den aufgeführten Befehlen und Beispielen orientiert sich an dem Inhaltsverzeichnis der 17. Auflage der ‘Angewandten Statistik’. Nähere Hinweise zu den verwendeten Formeln sowie Erklärungen zu den Beispielen sind in dem Buch (Ebook) nachzulesen!

Aktuelle R-Umgebung zu den Beispielen

The R Project for Statistical Computing

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.6.3  magrittr_1.5    tools_3.6.3     htmltools_0.4.0
##  [5] yaml_2.2.1      Rcpp_1.0.4      stringi_1.4.6   rmarkdown_2.1  
##  [9] knitr_1.28      stringr_1.4.0   xfun_0.12       digest_0.6.25  
## [13] rlang_0.4.5     evaluate_0.14

5.1 Zufallsvariable

Augenzahl beim Werfen eines Würfels:

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
x <- seq(1, 6, by=1); fx <- rep(1/6, 6)
plot(x, fx, type="h", ylim=c(0, 0.2), ylab="f(x)", 
                              xlab="Augenzahl", col="grey")
points(x, fx, pch=19, cex=1.1, col="black")
F <- cumsum(fx)
x1 <- c(0, x, 7)
Fx <- c(0, F, 1)
plot(x1, Fx, type="s", ylim=c(0, 1), ylab="F(x)", xlab="Augenzahl")

Augenzahl beim Werfen zweier Würfel:

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
x <- seq(2, 12, by=1)
fx <- c(1/36, 2/36, 3/36, 4/36, 5/36, 6/36, 5/36, 4/36, 3/36, 2/36, 1/36)
plot(x, fx, type="h", ylim=c(0, 0.2), ylab="f(x)", 
                                            xlab="Augenzahl", col="grey")
points(x, fx, pch=19, cex=1.1, col="black")
F <- cumsum(fx)
x1 <- c(0, x, 13)
Fx <- c(0, F, 1)
plot(x1, Fx, type="s", ylim=c(0,   1), ylab="F(x)", xlab="Augenzahl")

5.2 Maßzahlen zur Kennzeichnung einer Verteilung

5.2.3.1 Momente: Schiefe und Exzess

Funktionen skewness() und kurtosis() in library(e1071)

x <- c(2, 3, 4, 4, 4, 5, 5, 5, 5, 6, 8, 10, 20, 40)                                                     
m <- mean(x); s <- sd(x); n <- length(x)

(1/n)*sum((x-m)^3)/(s^3)                                 # Schiefe 
## [1] 2.198071

(1/n)*sum((x-m)^4)/(s^4) - 3                             # Wölbung
## [1] 3.89879

library(e1071)
skewness(x, type = 3)                                    # Schiefe    
## [1] 2.198071

kurtosis(x, type = 3)                                    # Wölbung     
## [1] 3.89879

5.2.3.2 Potenzmomente aus klassierten Häufigkeiten

momente <- function(x, f, d, b) {                        # x Daten (Vektor)
                                                         # f Klassenhäufigkeiten
                                                         # d Ursprung (Modalwert, Referenz)
                                                         # b Klassenbreite
    n <- sum(f);     z <- (x - d) / b
    m1 <- sum(f * z) / n;      m2 <- sum(f * z^2) / n
    m3 <- sum(f * z^3) / n;    m4 <- sum(f * z^4) / n
    s2 <- b^2 * (m2 - m1^2);   s3 <- (sqrt(s2))^3;     s4 <- s2^2
    mitt <- d + (b * m1)
    vari <- b^2 * (m2 - m1^2)
    schi <- (b^3 * (m3 - 3*m1*m2 + 2*m1^3)) / s3
    woel <- ((b^4 * (m4 - 4*m1*m3 + 6*m1^2*m2 - 3*m1^4)) / s4) - 3
    cat("\n","Momente aus klassierten Beobachtungen:","\n",
        "Mittelwert:",mitt,"\n",        "   Varianz:",vari,"\n",
        "   Schiefe:",schi,"\n",        "   Wölbung:",woel,"\n")
} 

x <- c(8.8, 9.3, 9.8, 10.3, 10.8, 11.3, 11.8)          
f <- c(  4,   8,  11,    7,    5,    3,    2)
d <- 9.8;   b <- 0.5

momente(x, f, d, b)                                                
## 
##  Momente aus klassierten Beobachtungen: 
##  Mittelwert: 10.025 
##     Varianz: 0.636875 
##     Schiefe: 0.4598458 
##     Wölbung: -0.4809175

Unterschiedliche Verteilungsformen:

b      <- 10                                             # Klassenbreite
d      <- 30                                             # Ursprung
x      <- c(15, 25, 35, 45, 55, 65, 75, 85)              # Klassenmitten

f1     <- c( 5, 10, 15, 30, 30, 15, 10,  5)              # Häufigkeiten
momente(x, f1, d, b); barplot(f1, width=b, space=0, names.arg=x)
## 
##  Momente aus klassierten Beobachtungen: 
##  Mittelwert: 50 
##     Varianz: 275 
##     Schiefe: 0 
##     Wölbung: -0.3140496


f2     <- c(20, 30, 30, 20, 10,  5,  5,  0)              # Häufigkeiten
momente(x, f2, d, b); barplot(f2, width=b, space=0, names.arg=x)
## 
##  Momente aus klassierten Beobachtungen: 
##  Mittelwert: 35.41667 
##     Varianz: 245.6597 
##     Schiefe: 0.7101991 
##     Wölbung: -0.02114747


f3     <- c( 0,  5,  5, 10, 20, 30, 30, 20)              # Häufigkeiten
momente(x, f3, d, b); barplot(f3, width=b, space=0, names.arg=x)
## 
##  Momente aus klassierten Beobachtungen: 
##  Mittelwert: 64.58333 
##     Varianz: 245.6597 
##     Schiefe: -0.7101991 
##     Wölbung: -0.02114747


f4     <- c( 5, 15, 30, 10, 10, 30, 15,  5)              # Häufigkeiten
momente(x, f4, d, b); barplot(f4, width=b, space=0, names.arg=x)
## 
##  Momente aus klassierten Beobachtungen: 
##  Mittelwert: 50 
##     Varianz: 375 
##     Schiefe: 0 
##     Wölbung: -1.235556

5.2.3.3 Quantilsmaße zu Schiefe und Exzess

                                                         # Körpergröße von 70 Studenten
y <- c(63, 63, 64, 64, rep(65,4), rep(66,5), rep(67,4), 
       rep(68,6), rep(69,5), rep(70,8), rep(71,7), rep(72,7), 
       rep(73,10), rep(74,5), rep(75,3), rep(76,2))

par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=11)
hist(y, breaks=c(seq(62, 76, by=2)), freq=TRUE, col="grey",
     ylim=c(0, 15), xlim=c(62, 76), xlab=" ", main=" ", ylab="Häufigkeit")


mean(y)                                                  # Mittelwert
## [1] 70.04286
var(y)                                                   # empirische Varianz
## [1] 11.11408
skewness(y)                                              # empirisches 3tes Moment
## [1] -0.2843902
kurtosis(y)                                              # empirisches 4tes Moment
## [1] -0.8728042

                                                         # Quartile
Q    <- quantile(y, probs=seq(0, 1, 0.25), names=TRUE, type=7); Q    
##   0%  25%  50%  75% 100% 
##   63   68   70   73   76

Q1   <- as.numeric(Q[2]); Q2 <- as.numeric(Q[3]); Q3 <- as.numeric(Q[4])
skew <- (Q3 + Q1 - 2*Q2)/(Q3-Q1);  skew
## [1] 0.2
                                                         # Oktile  
A    <- quantile(y, probs=seq(0, 1, 0.125), names=TRUE, type=7); A  
##    0% 12.5%   25% 37.5%   50% 62.5%   75% 87.5%  100% 
##    63    66    68    69    70    72    73    74    76
A7 <- as.numeric(A[8]); A6 <- as.numeric(A[7]); A5 <- as.numeric(A[6])
A3 <- as.numeric(A[4]); A2 <- as.numeric(A[3]); A1 <- as.numeric(A[2])
kurt <- ((A7 - A5) + (A3 - A1))/(A6-A2); kurt
## [1] 1

Tukey’s Fünfer-Regel:

                                                         # Körpergröße von 70 Studenten
y <- c(63, 63, 64, 64, rep(65,4), rep(66,5), rep(67,4), 
       rep(68,6), rep(69,5), rep(70,8), rep(71,7), rep(72,7), 
       rep(73,10), rep(74,5), rep(75,3), rep(76,2))
fivenum(y)                                               # Tukeys five numbers
## [1] 63 68 70 73 76

5.3 Diskrete Verteilungen

5.3.2 Gleichverteilung

Funktionen zu diskreten Verteilungen in library(e1071)

library(e1071)
ddiscrete(1:10, rep(0.1, 10))                            # Dichtefunktion             
##  [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
pdiscrete(1:10, rep(0.1, 10))                            # Verteilungsfunktion        
##  [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
qdiscrete(c(0.25, 0.5, 0.75), rep(0.1, 10))              # Quantilfunktion (Quartile) 
## [1] 3 5 8
rdiscrete(20, 1:10)                                      # Zufallszahlen              
##  [1]  7  7  9  8  3  5  7  7 10  7  3  4  9 10  3 10  9  5  3  7

x  <- 1:10

fx <- ddiscrete(1:10, rep(0.1, 10))                      # Wahrscheinlichkeitsfunktion

Fx <- pdiscrete(1:10, rep(0.1, 10))                      # Verteilungsfunktion
x1 <- c(0, x, 11, 12)
Fx <- c(0, Fx, 1, 1)

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=10)
plot(x, fx, type="h", ylim=c(0, 0.12), xlim=c(0,12), ylab="f(x)", 
            xlab=" ", lty=2, col="grey")
points(x, fx, pch=19, cex=1.1, col="black")
plot(x1, Fx, type="s", ylim=c(0,1), xlim=c(0,12), ylab="F(x)", xlab=" ")

5.3.3 Binomialverteilung

Bernoulli Versuch:

n <- 100
x <- sample(c(-1,1), n, replace=T, prob=c(0.6,0.4))
                                                          
par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12)
plot(x, type="h", axes=FALSE, xlab="Bernoulli Versuche", ylab="Ausgang")
axis(1, at = NULL, labels = TRUE, tick = TRUE)
text(0, -0.5, "A");  text(0, +0.5, "B")

Binomialverteilung:

n  <- 4; p  <- 1/6;  x  <- 0:n 

fx <- dbinom(x, n, p)                                     # Wahrscheinlichkeitsfunktion

Fx <- pbinom(x, n, p); x1 <- c(x, n+1, n+2)               # Verteilungsfunktion
Fx <- c(Fx, 1, 1)

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=10)
plot(x, fx, type="h", ylim=c(0, 0.5), xlim=c(0,n), ylab="f(x)", 
            xlab=" ", lty=2, col="grey")
points(x, fx, pch=19, cex=1.1, col="black")
plot(x1, Fx, type="s", ylim=c(0,1), xlim=c(0,n+2), ylab="F(x)", xlab=" ")
lines(c(0,0),c(0,Fx[1]))

dbinom(3, 3, 1/2)
## [1] 0.125
dbinom(2, 3, 1/2)
## [1] 0.375

Tabelle zu Binomialwahrscheinlichkeiten:

prob <- c(0.01, 0.05, 0.10, 0.20, 0.25, 0.30, 0.40, 0.50)
tab  <- matrix(data = NA, nrow = 18, ncol = 8, byrow = FALSE, dimnames = NULL)
row  <- 1
for (n in 2:5) { 
  for (i in 0:n) {
    tab[row, ] <- round(dbinom(i, n, prob), 4)
    row <- row +1  }  }
c1 <- c(rep(2,3), rep(3,4), rep(4,5), rep(5,6))
c2 <- c(0:2, 0:3, 0:4, 0:5)
tab <- cbind(c1, c2, tab) 
colnames(tab) <- c("n", "x", "p=0.01", "p=0.05", "p=0.10",  "p=0.20", "p=0.25", "p=0.30", "p=0.40", "p=0.50")
as.data.frame(tab)  

Beispiele zur Binomialverteilung:

dbinom(0, 4, 0.2)                                           # Ausschussware
## [1] 0.4096
dbinom(1, 4, 0.2)
## [1] 0.4096
dbinom(2, 4, 0.2)
## [1] 0.1536

dbinom(0:4, 4, 0.2)
## [1] 0.4096 0.4096 0.1536 0.0256 0.0016
pbinom(2, 4, 0.2)
## [1] 0.9728
                                         
1 - pbinom(0, 6, 1/6, lower.tail=TRUE)                      # Chevalier de Mere
## [1] 0.665102
pbinom(1, 12, 1/6, lower.tail=FALSE)
## [1] 0.6186674
                                          
pbinom(18, 120, 1/6)                                        # Würfel
## [1] 0.3657008
                                      
round(200*dbinom(0:4, 4, 0.465),  2)                        # Mäuse
## [1] 16.38 56.96 74.27 43.03  9.35
                                            
dbinom(1, 2, 0.8)                                           # Behandlungserfolge
## [1] 0.32
dbinom(1, 5, 0.8)
## [1] 0.0064
dbinom(5, 5, 0.8)
## [1] 0.32768
                                            
dbinom(2, 5, 0.5)                                           # fünf Kinder
## [1] 0.3125
dbinom(5, 5, 0.5)
## [1] 0.03125

5.3.3.3 Approximation durch die Standardnormalverteilung

Tabelle zum Vergleich unterschiedlicher Tranformationen:

p <- c(0.05, 0.5); q <- 1-p

n <- c(100,   10)
X <- c(10,     8)
                                                        
P0 <- 1 - pbinom(X, n, p, lower.tail = TRUE)                # exakt binomial  

m <- n*p; s <- sqrt(n*p*(1-p))
Z1 <- (X + 0.5 - m) / s                                     # Normalverteilung
P1 <- 1 - pnorm(Z1)

Z2 <- (asin(sqrt(X/n)) - asin(sqrt(p))) / sqrt(1/(4*n))     # ArcSin-Transf. 
P2 <- 1 - pnorm(Z2)

Z3 <- 2*(asin(sqrt((X+1)/(n+1))) - asin(sqrt(p))) * sqrt(n+0.5)  # Freeman
P3 <- 1 - pnorm(Z3)

Z4 <- (sqrt(4*q*(X+1)) - sqrt(4*p*(n-X)))                   # Mosteller-Tukey 
P4 <- 1 - pnorm(Z4)

Z5 <- (sqrt(q*(4*X+3)) - sqrt(p*(4*n-4*X-1)))               # Freeman-Tukey   
P5 <- 1 - pnorm(Z5)

tab <- round(cbind(P0, P1, P2, P3, P4, P5), 4)        
c1  <- c(0.05, 0.50); c2  <- c(100, 10); c3 <- c(10, 8)
tab <- cbind(c1, c2, c3, tab)
colnames(tab) <- c("p", "n", "X", "exakt ", "normal",  "arcsin", "Freeman", "Mosteller/Tukey", "Freeman/Tukey")
as.data.frame(tab)

Tabelle zu Winkeltransformationen:

p.1 <- seq(0.0, 0.9, by=0.1)                        
p.2 <- seq(0, 0.09, by=0.01)
p   <- matrix(rep(NA, 100), nrow=10, ncol=10)
for (i in 1:10) {for (k in 1:10) {p[i,k] <- p.1[i] + p.2[k]}}
z   <- round(asin(sqrt(p))*(180/pi), 3)
tab <- cbind(p.1, z)
colnames(tab) <- c("p", seq(0.00, 0.09, 0.01))
as.data.frame(tab)

5.3.4 Multinomialverteilung

Funktion scatterplot3d() in library(scatterplot3d)

n <- 4; p=c(0.4, 0.3, 0.3)
M <- matrix(c(0,0,4, 0,1,3, 0,2,2, 0,3,1, 0,4,0, 1,0,3, 1,1,2, 1,2,1, 1,3,0,
              2,0,2, 2,1,1, 2,2,0, 3,0,1, 3,1,0, 4,0,0),
              byrow=FALSE, nrow=3);  
P <- rep(NA, 15); for (i in 1:15) P[i] <- dmultinom(M[,i], prob=p)
round(P, 3)
##  [1] 0.008 0.032 0.049 0.032 0.008 0.043 0.130 0.130 0.043 0.086 0.173 0.086
## [13] 0.077 0.077 0.026
                                                          
library(scatterplot3d)                                   
par(mfrow=c(1,1), lwd=2, font.axis=2, bty="n", ps=12)
z <- P; x <- M[1,]; y <- M[2,]
scatterplot3d(x, y, z, col.axis="black", type="h", highlight.3d=TRUE,
      box=FALSE, col.grid="grey", pch=20, cex.symbols=2,
      xlab="x (Anzahl für A)", ylab="y (Anzahl für B)", zlab="P(X=x, Y=y)")

Beispiele zur Multinomialverteilung:

dmultinom(c(3,2,1), prob=c(0.50, 0.30, 0.20))             # Beispiel Perlen 
## [1] 0.135

dmultinom(c(1,1,1,3,3,3), prob=rep(1/6, 6))               # Beispiel Würfel
## [1] 0.001018751
                                                           
dmultinom(c(8,1,1), prob=c(1/3, 1/3, 1/3))                # Kandidatenwahl 
## [1] 0.001524158
dmultinom(c(3,3,4), prob=c(1/3, 1/3, 1/3))
## [1] 0.07112737

5.3.5 Poisson-Verteilung

x  <- 0:12
f1 <- dpois(x, 1); f2 <- dpois(x, 2); f3 <- dpois(x, 6)

par(mfrow=c(1,3), lwd=2, font.axis=2, bty="n", ps=14)
plot(x, f1, type="h", ylim=c(0, 0.4), xlim=c(0,12),
     ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(x, f1, pch=19, cex=1.2, col="black")
text(6, 0.18, expression(lambda == 1), cex=1.5) 
plot(x, f2, type="h", ylim=c(0, 0.3), xlim=c(0,12),
     ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(x, f2, pch=19, cex=1.2, col="black") 
text(6, 0.28, expression(lambda == 2), cex=1.5) 
plot(x, f3, type="h", ylim=c(0, 0.2), xlim=c(0,12),
     ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(x, f3, pch=19, cex=1.2, col="black") 
text(6, 0.18, expression(lambda == 6), cex=1.5) 

Tabelle zur Poissonverteilung:

lambda <- c(0.2, 0.5, 0.8, 1, 3, 5, 8, 12, 20)
tab <- matrix(data = NA, nrow = 30, ncol = 9, 
              byrow = FALSE, dimnames = NULL)

for (i in 0:29) {tab[i+1, ] <- round(dpois(i, lambda), 4) }

tab <- cbind(0:29, tab)
colnames(tab) <- c("x","lambda=0.2", "lambda=0.5", "lambda=0.8", 
        "lambda=1", "lambda=3", "lambda=5", "lambda=8", "lambda=12", "lambda=20")
as.data.frame(tab[,1:5])  

Beispiele zur Poissonverteilung:

dpois(0:3, 2.7397)                                           # Geburtstagsproblem
## [1] 0.06458972 0.17695646 0.24240380 0.22137123

                                               
dpois(3, 2)                                                  # Serumproben       
## [1] 0.180447
1-ppois(2, 2,)
## [1] 0.3233236

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
plot(0:10, dpois(0:10, 2), type="h", ylim=c(0, 0.3), xlim=c(0,10),
     ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(0:10, dpois(0:10, 2), pch=19, cex=1.0, col="black") 
plot(0:7, ppois(0:7, 2), type="s", ylim=c(0,1), 
    xlim=c(0,10), ylab="F(x)", xlab=" ") 
lines(c(0,0), c(0,dpois(0, 2)))

5.3.5.1 Dispersionsindex

DI_test <- function(xi, fi) {
  if (length(xi)!=length(fi)) stop("x unf f ungleiche Längen!!!") 
  n    <- sum(fi);    mw   <- sum(xi*fi)/n
  var  <- (sum(xi^2*fi) - (sum(xi*fi))^2/n) / (n-1)
  DI   <- mw/var
  stat <- sum(fi*(xi-mw)^2)/mw
  pvalue <- pchisq(stat, df=n-1, lower.tail=F)
  cat("\n","Mittelwert:",round(mw, 4),"\n","   Varianz:",round(var,4),
      "\n","        DI:",round(DI, 4),
      "\n","  Statstik:",stat,"(P-Wert:", pvalue,")","\n")     }

                                                             # Tod durch Pferdehufschlag 
x <- c(  0,  1,  2, 3, 4, 5)                  
f <- c(109, 65, 22, 3, 1, 0)
DI_test(x, f)
## 
##  Mittelwert: 0.61 
##     Varianz: 0.611 
##          DI: 0.9984 
##    Statstik: 199.3115 (P-Wert: 0.4804495 )

lambda <- 0.61
n      <- 200                                
round(dpois(0:5, lambda) * n, 1)
## [1] 108.7  66.3  20.2   4.1   0.6   0.1

5.3.5.2 Approximation durch die Standardnormalverteilung

k <- c(3, 4); lambda <- c(9, 10)                             # Approximationen
z1 <- (k  - lambda) / sqrt(lambda); z
##  [1] 0.0081 0.0324 0.0486 0.0324 0.0081 0.0432 0.1296 0.1296 0.0432 0.0864
## [11] 0.1728 0.0864 0.0768 0.0768 0.0256
ppois(k, lambda)
## [1] 0.02122649 0.02925269
pnorm(z1)
## [1] 0.02275013 0.02888979

5.3.6 Negative Binomialverteilung

choose(7+3-1, 7) * 0.2^3 * 0.8^7
## [1] 0.06039798
dnbinom(7, 3, 0.2)
## [1] 0.06039798

p <- rep(NA, 8)
for (i in 0:7) p[i+1] <- choose(i+3-1, i) * 0.2^3 * 0.8^i; sum(p)
## [1] 0.3222005
pnbinom(7, 3, 0.2)
## [1] 0.3222005
par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
plot(0:40, dnbinom(0:40, 3, 0.2), type="h", ylim=c(0, 0.08), 
    xlim=c(0,40), ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(0:40, dnbinom(0:40, 3, 0.2), pch=19, cex=0.5, col="black")
plot(0:40, pnbinom(0:40, 3, 0.2), type="s", ylim=c(0,1), 
    xlim=c(0,40), ylab="F(x)", xlab=" ")
lines(c(0,0), c(0,pnbinom(0, 3,0.2)))

Beispiele zur negativen Binomialverteilung:

p <- 0.4; x <- 10; k <- 5                                   # Werfen mit Steinen
choose(x-1, k-1)*p^k*(1-p)^(x-k)
## [1] 0.1003291

dnbinom(x-k, k, 0.4)
## [1] 0.1003291
                                                        
k   <- c(  0,   1,  2,  3, 4, 5)                            # Unfälle
obs <- c(447, 132, 42, 21, 3, 2); n <- sum(obs)
m   <- sum(obs*k)/n; round(m, 2)                            # Mittelwert (Erwartungswert)   
## [1] 0.47
round(dpois(k, m) * n, 0)                                   # Poisson- Verteilung           
## [1] 406 189  44   7   1   0
v   <- sum((obs*(k - m)^2))/(n-1); v                        # (emp.) Varianz                
## [1] 0.6919002
p   <- m / v;  r <- m*p/(1-p)                               # Modellparameter               
round(dnbinom(k, r, p) * n, 0)                              # negative Binomialverteilung   
## [1] 443 139  44  14   5   2

                                                            #  Zecken
beob <- c(rep(0,7), rep(1,9), rep(2,8), rep(3,13), rep(4,8), rep(5,5), 
          rep(6,4), rep(7,3), rep(8,0), rep(9,1), 10, 10); n  <- sum(beob) 
r.hat <- mean(beob)^2/(var(beob)-mean(beob)); r.hat
## [1] 3.956746
p.hat <- r.hat/(mean(beob)+r.hat); p.hat
## [1] 0.5490336
round(dnbinom(0:11, 3.96, 0.55)*60, 0)                                                 
##  [1]  6 10 11 10  8  6  4  2  1  1  1  0

                                                            # Markenartikel
beob <- c(rep(0,39), rep(1,14), rep(2,10), rep(3,6), rep(4,4), rep(5,4),
          rep(6,3), rep(7,3), rep(8,2), rep(9,2), rep(10, 13))
n  <- sum(beob); m  <- mean(beob); m 
## [1] 2.91
v  <- var(beob)
r  <- m^2 / (v - m); r 
## [1] 0.8536217

m  <- 3.4;  s  <- 0.5;  p  <- s/(s+m)                    
n=100; x  <- 0:10
round(dnbinom(x, s, p)*n, 0)
##  [1] 36 16 10  7  6  4  4  3  2  2  2
round(dpois(x, m)*n, 0)
##  [1]  3 11 19 22 19 13  7  3  1  1  0

5.3.7 Geometrische Verteilung

p <- 1/6; n <- 0:20
f <- dgeom(n, p)                                             # Wahrscheinlichkeitsfunktion
F <- pgeom(n, p)                                             # Verteilungsfunktion
                                                        
par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
plot(n, f, type="h", ylim=c(0, 0.2), 
        xlim=c(0,20), ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(n, f, pch=19, cex=0.5, col="black")
plot(n, F, type="s", ylim=c(0,1), xlim=c(0,20), ylab="F(x)", xlab=" ")
lines(c(0,0), c(0,p))

5.3.8 Hypergeometrische Verteilung

dhyper(2, 5, 10, 5)                                          # Beispiel Urnenmodell
## [1] 0.3996004
                                                      
dhyper(1:5, 6, 4, 5)                                         # Beispiel Studenten
## [1] 0.02380952 0.23809524 0.47619048 0.23809524 0.02380952
phyper(1:5, 6, 4, 5)
## [1] 0.02380952 0.26190476 0.73809524 0.97619048 1.00000000
                                                        
dhyper(4, 6, 43, 6)                                          # Beispiel Lotto
## [1] 0.0009686197
                                                           
x <- 0:6
f <- dhyper(0:6, 6, 43, 6)                                   # Wahrscheinlichkeitsfunktion
F <- phyper(0:6, 6, 43, 6)                                   # Verteilungsfunktion

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11)
plot(x, f, type="h", ylim=c(0, 0.5), xlim=c(0,7),
             ylab="f(x)", xlab=" ", lty=2, col="grey", lwd=1)
points(x, f, pch=19, cex=1.0, col="black") 
plot(x, F, type="s", ylim=c(0,1), xlim=c(0,7), ylab="F(x)", xlab=" ") 
lines(c(0,0), c(0,dhyper(0, 6, 43, 6)))