1 R-Pakete

library(tidyverse)  # data wrangling
theme_set(theme_minimal())  # Stylesheet für ggplot2

2 Hintergrund

Will man eine Verteilung untersuchen, sind Verteilungsfunktion $F$ und Quantilsfunktion $F^{- 1}$ wichtige Größen. Nicht nur für theoretische, sondern auch für empirische Verteilungen kann man diese Funktionen anwenden.

Allerdings ist die Umsetzung in R vielleicht nicht ganz klar. Daher soll dieser Post zeigen, wie man eine empirische Verteilungfunktion und Quantilsfunktion in R berechnet.

3 Verteilungsfunktion der Normalverteilung

Für eine Normalverteilung kann man sich Quantile und Verteilungsfunktion einfach in R ausgeben lassen.

Wert der Verteilungsfunktion bei $x = 0$ :

pnorm(q = 0)
#> [1] 0.5

Ausgehend von einer Standardnormalverteilung sagt uns R, dass der Wert der Verteilungsfunktion $F (0) = 0.5$ beträgt.

Quantile, z.B. $F^{- 1} (0.5)$ :

qnorm(p = .5)
#> [1] 0

Passt!

Das Plotten geht analog:

ggplot(tibble(x = c(-3, 3)), aes(x = x)) +
  stat_function(fun = pnorm) +
  labs(title = "Verteilungsfunktion F der Normalverteilung")

4 Empirische Verteilungsfunktion

Erzeugen wir uns (standard-)normalverteilte Daten:

d_normal <-
  tibble(x = rnorm(n = 1e04))

4.1 Tidyverse

Die empirische Verteilungsfunktion für $x = 0$ kann man, so mit dem Tidyverse bekommen

4.1.1 Tidyverse 1

Dann zählen wir den Anteil der Beobachtungen (Zeilen), die dem gesuchten Wert der Verteilungsfunktion entsprechen, z.B. $x = 0$ :

d_normal %>% 
  count(x <= 0)
#> # A tibble: 2 × 2
#>   `x <= 0`     n
#>   <lgl>    <int>
#> 1 FALSE     4951
#> 2 TRUE      5049

Also etwa 50%.

4.1.2 Tidyverse 2

So geht es auch, recht umständlich allerdings:

Zunächst erzeugen wir die kumulierte Verteilung.

d_normal <-
  d_normal %>% 
  mutate(x_cume_dist = cume_dist(x))

Dann lesen wir den Wert von x_cume_dist an der gewünschten Stelle von x aus:

d_normal %>% 
  filter(x <= 0) %>% 
  summarise(max(x_cume_dist))
#> # A tibble: 1 × 1
#>   `max(x_cume_dist)`
#>                <dbl>
#> 1              0.505

4.1.3 Plotten der ECDF

Und so kann man die empirische Verteilungsfunktion plotten, auch genannt empirical cumulative density function (ecdf):

d_normal %>%
  ggplot(aes(x = x)) +
  stat_ecdf(geom = "step")

4.1.4 Quantile

Zur Erinnerung: Das $p$ -Quantil ist der Wert, der von $p \cdot 100$ Prozent der Beobachtungen der Verteilung nicht überschritten wird.

Anschaulich: Wir schneiden den Anteil $p$ links von der Verteilung ab, und fragen uns, bei welchem Wert $x$ wir die Schere ansetzen müssen.

Quantil für $p = .5$ :

d_normal %>% 
  filter(percent_rank(x) <= .5) %>% 
  summarise(max(x))
#> # A tibble: 1 × 1
#>   `max(x)`
#>      <dbl>
#> 1  -0.0110

Also etwa Null, passt.

4.2 Base R

4.2.1 Quantile

Mit Base R geht es schön einfach:

quantile(d_normal$x, prob = .5)
#>         50% 
#> -0.01070173

Oder, etwas mehr Tidyverse-Stil:

d_normal %>% 
  summarise(x_50perc = quantile(x, prob = .5))
#> # A tibble: 1 × 1
#>   x_50perc
#>      <dbl>
#> 1  -0.0107

4.2.2 ECDF

mit ecdf(x) erzeugt man eine Funktion und zwar eine, die die kumulierten Anteile für x wiedergibt:

F <- ecdf(d_normal$x)

Wenn einmal definiert, können wir die Funktion bequem befragen.

Hey F, was ist dein Wert bei $x = 0$ ?

F(0)
#> [1] 0.5049

50%! Passt

4.2.3 Plot

ecdf() hat auch eine Plot-Methode:

plot(ecdf(d_normal$x))

4.3 Mosaic

library(mosaic)

Mit ${{m o s a i c}}$ ` geht das auch schön einfach:

4.3.1 ECDF

Was ist der kumulierte Anteil für x=0:

pdata( ~ x, q = 0, data = d_normal)
#> [1] 0.5049

4.3.2 Quantile

Was ist das Quantil für p = .5:

qdata(~ x, p = 1/2, data = d_normal)
#>         50% 
#> -0.01070173

Etwa Null, passt!

5 Reproducibility

#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.1.3 (2022-03-10)
#>  os       macOS Big Sur/Monterey 10.16
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Europe/Berlin               
#>  date     2022-05-02                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  package     * version date       lib source                            
#>  assertthat    0.2.1   2019-03-21 [1] CRAN (R 4.1.0)                    
#>  backports     1.4.1   2021-12-13 [1] CRAN (R 4.1.0)                    
#>  blogdown      1.8     2022-02-16 [2] CRAN (R 4.1.2)                    
#>  bookdown      0.24.2  2021-10-15 [1] Github (rstudio/bookdown@ba51c26) 
#>  brio          1.1.3   2021-11-30 [1] CRAN (R 4.1.0)                    
#>  broom         0.7.12  2022-01-28 [1] CRAN (R 4.1.2)                    
#>  bslib         0.3.1   2021-10-06 [1] CRAN (R 4.1.0)                    
#>  cachem        1.0.6   2021-08-19 [1] CRAN (R 4.1.0)                    
#>  callr         3.7.0   2021-04-20 [1] CRAN (R 4.1.0)                    
#>  cellranger    1.1.0   2016-07-27 [1] CRAN (R 4.1.0)                    
#>  cli           3.2.0   2022-02-14 [1] CRAN (R 4.1.2)                    
#>  codetools     0.2-18  2020-11-04 [2] CRAN (R 4.1.3)                    
#>  colorout    * 1.2-2   2022-01-04 [1] Github (jalvesaq/colorout@79931fd)
#>  colorspace    2.0-3   2022-02-21 [1] CRAN (R 4.1.2)                    
#>  crayon        1.5.1   2022-03-26 [1] CRAN (R 4.1.2)                    
#>  DBI           1.1.2   2021-12-20 [1] CRAN (R 4.1.0)                    
#>  dbplyr        2.1.1   2021-04-06 [1] CRAN (R 4.1.0)                    
#>  desc          1.4.0   2021-09-28 [1] CRAN (R 4.1.0)                    
#>  devtools      2.4.3   2021-11-30 [1] CRAN (R 4.1.0)                    
#>  digest        0.6.29  2021-12-01 [1] CRAN (R 4.1.0)                    
#>  dplyr       * 1.0.8   2022-02-08 [1] CRAN (R 4.1.2)                    
#>  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.1.0)                    
#>  evaluate      0.14    2019-05-28 [1] CRAN (R 4.1.0)                    
#>  fansi         1.0.3   2022-03-24 [1] CRAN (R 4.1.2)                    
#>  fastmap       1.1.0   2021-01-25 [2] CRAN (R 4.1.0)                    
#>  forcats     * 0.5.1   2021-01-27 [1] CRAN (R 4.1.0)                    
#>  fs            1.5.2   2021-12-08 [1] CRAN (R 4.1.0)                    
#>  generics      0.1.2   2022-01-31 [1] CRAN (R 4.1.2)                    
#>  ggplot2     * 3.3.5   2021-06-25 [2] CRAN (R 4.1.0)                    
#>  glue          1.6.2   2022-02-24 [1] CRAN (R 4.1.2)                    
#>  gtable        0.3.0   2019-03-25 [1] CRAN (R 4.1.0)                    
#>  haven         2.4.3   2021-08-04 [1] CRAN (R 4.1.0)                    
#>  hms           1.1.1   2021-09-26 [1] CRAN (R 4.1.0)                    
#>  htmltools     0.5.2   2021-08-25 [1] CRAN (R 4.1.0)                    
#>  httr          1.4.2   2020-07-20 [1] CRAN (R 4.1.0)                    
#>  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.1.0)                    
#>  jsonlite      1.7.3   2022-01-17 [1] CRAN (R 4.1.2)                    
#>  knitr         1.37    2021-12-16 [1] CRAN (R 4.1.0)                    
#>  lifecycle     1.0.1   2021-09-24 [1] CRAN (R 4.1.0)                    
#>  lubridate     1.8.0   2021-10-07 [1] CRAN (R 4.1.0)                    
#>  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.1.2)                    
#>  memoise       2.0.0   2021-01-26 [2] CRAN (R 4.1.0)                    
#>  modelr        0.1.8   2020-05-19 [1] CRAN (R 4.1.0)                    
#>  munsell       0.5.0   2018-06-12 [1] CRAN (R 4.1.0)                    
#>  pillar        1.7.0   2022-02-01 [1] CRAN (R 4.1.2)                    
#>  pkgbuild      1.2.0   2020-12-15 [2] CRAN (R 4.1.0)                    
#>  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.1.0)                    
#>  pkgload       1.2.4   2021-11-30 [1] CRAN (R 4.1.0)                    
#>  prettyunits   1.1.1   2020-01-24 [1] CRAN (R 4.1.0)                    
#>  processx      3.5.2   2021-04-30 [1] CRAN (R 4.1.0)                    
#>  ps            1.6.0   2021-02-28 [1] CRAN (R 4.1.0)                    
#>  purrr       * 0.3.4   2020-04-17 [1] CRAN (R 4.1.0)                    
#>  R6            2.5.1   2021-08-19 [1] CRAN (R 4.1.0)                    
#>  Rcpp          1.0.8.3 2022-03-17 [1] CRAN (R 4.1.2)                    
#>  readr       * 2.1.2   2022-01-30 [1] CRAN (R 4.1.2)                    
#>  readxl        1.3.1   2019-03-13 [1] CRAN (R 4.1.0)                    
#>  remotes       2.4.0   2021-06-02 [2] CRAN (R 4.1.0)                    
#>  reprex        2.0.1   2021-08-05 [1] CRAN (R 4.1.0)                    
#>  rlang         1.0.2   2022-03-04 [1] CRAN (R 4.1.2)                    
#>  rmarkdown     2.11    2021-09-14 [1] CRAN (R 4.1.0)                    
#>  rprojroot     2.0.2   2020-11-15 [2] CRAN (R 4.1.0)                    
#>  rstudioapi    0.13    2020-11-12 [1] CRAN (R 4.1.0)                    
#>  rvest         1.0.2   2021-10-16 [1] CRAN (R 4.1.0)                    
#>  sass          0.4.0   2021-05-12 [1] CRAN (R 4.1.0)                    
#>  scales        1.2.0   2022-04-13 [1] CRAN (R 4.1.3)                    
#>  sessioninfo   1.1.1   2018-11-05 [2] CRAN (R 4.1.0)                    
#>  stringi       1.7.6   2021-11-29 [1] CRAN (R 4.1.0)                    
#>  stringr     * 1.4.0   2019-02-10 [1] CRAN (R 4.1.0)                    
#>  testthat      3.1.2   2022-01-20 [1] CRAN (R 4.1.2)                    
#>  tibble      * 3.1.6   2021-11-07 [1] CRAN (R 4.1.0)                    
#>  tidyr       * 1.2.0   2022-02-01 [1] CRAN (R 4.1.2)                    
#>  tidyselect    1.1.2   2022-02-21 [1] CRAN (R 4.1.2)                    
#>  tidyverse   * 1.3.1   2021-04-15 [1] CRAN (R 4.1.0)                    
#>  tzdb          0.1.2   2021-07-20 [2] CRAN (R 4.1.0)                    
#>  usethis       2.0.1   2021-02-10 [2] CRAN (R 4.1.0)                    
#>  utf8          1.2.2   2021-07-24 [1] CRAN (R 4.1.0)                    
#>  vctrs         0.4.0   2022-03-30 [1] CRAN (R 4.1.2)                    
#>  withr         2.5.0   2022-03-03 [1] CRAN (R 4.1.2)                    
#>  xfun          0.29    2021-12-14 [1] CRAN (R 4.1.0)                    
#>  xml2          1.3.3   2021-11-30 [1] CRAN (R 4.1.0)                    
#>  yaml          2.2.2   2022-01-25 [1] CRAN (R 4.1.2)                    
#> 
#> [1] /Users/sebastiansaueruser/Library/R/x86_64/4.1/library
#> [2] /Library/Frameworks/R.framework/Versions/4.1/Resources/library