knitr::opts_chunk$set(echo = TRUE,fig.width=6,fig.height=4)
    library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
    library(flextable)
## 
## Attaching package: 'flextable'
## 
## The following object is masked from 'package:purrr':
## 
##     compose
    library(ggplot2)
    library(dplyr)
    theme_set(theme_minimal()) # GGplot tema
    options(digits=3)
    pr <- function(data) {flextable(data)}

Anvisningar

I markdownfilen nedan har vi genomfört några exempel på hanteringen av diskreta bivariata fördelningar som tabeller.

Men det finns också stycken markerade med Uppgift med uppgifter att lösa med hjälp av R. Lösningen kan oftast hämtas mer eller mindre direkt ur exemplet som föregår.

Skriv in era lösningar i ett markdown-dokument och typsätt till en fil. Ni kan som mall använda detta dokument där avsnitten från Anvisningar till och med Betingade fördelningar har tagits bort.

En bivariat fördelning

Här definierar vi en tabell med kolumnerna X och Y samt en kolumn p(x,y) för den simultana massfunktionen \(p(x,y)=P(Y=y,X=x)\).

Med crossing() skapas en tabell vars rader svarar mot alla \((x,y)\)-värden produktmängden \(S_X \times S_Y\). Observera hur vi använder “backticks” (grav accent) för namnet p(x,y) som innehåller parenteser och komman.

    bivar1 <-
      crossing(X=1:3,Y=2:5) %>%
      mutate(`p(x,y)`=ifelse(X==Y,X,X+Y)) %>%
      mutate(`p(x,y)`=`p(x,y)`/sum(`p(x,y)`)) # Normalisera
    # Skriv ut
    pr(bivar1)

Ofta ges beskrivs en bivariat fördelning i en tabell med x- och y-värden givna i kolumner och rader istället. En så kallad “korstabell”. Man kan konvertera mellan representationerna med hjälp av pivot_wider() och pivot_longer().

Funktionen pivot_wider skapar en tabell där kolumnnamn hämtas från en kolumn och där värdena i tabellen från en annan.

    bivar1cross <-
      pivot_wider(bivar1,names_from = Y, values_from = `p(x,y)`)
    pr(bivar1cross)

Visualisering med ggplot

Att visualisera en multivariat fördelning diskret fördelning kan vara ganska komplicerat. Ett alternativ är att använda geom_tile() eller geom_raster().

    ggplot(bivar1, aes(x=X, y=Y, fill=`p(x,y)`)) +
      geom_raster() +
      geom_label(aes(label=round(`p(x,y)`,2)),color="white")

Beräkning av sannolikheter för händelser

Man kan beräkna sannolikheten för en händelse med hjälp av filter som väljer ut rader i tabellen och sedan summerar man med sum. Här beräknas sannolikheten för händelsen \(\{X+2 \ge Y\}\)

    filter(bivar1, X+2 >= Y) %>%
      summarise(`P(X+2>=Y)` = sum(`p(x,y)`))
## # A tibble: 1 × 1
##   `P(X+2>=Y)`
##         <dbl>
## 1       0.705

Man kan också notera att \(P((X,Y)\in A)\) kan skrivas som väntevärdet \[ E( \mathsf 1_{A} ) = \sum_{x,y} \mathsf 1_{A}(x,y) \cdot p(x,y) \] där \(\mathsf 1_A\) anger indikatorfunktionen för händelsen \(A\).

Här är motsvarande R-kod.

    bivar1 %>%
      mutate(A=ifelse((X+2) > Y,1,0)) %>%
      summarise(`P(A)` = sum(`p(x,y)` * A))
## # A tibble: 1 × 1
##   `P(A)`
##    <dbl>
## 1  0.410

Man kan åskådliggöra en händelse i XY-planet med geom_raster() och fill. Här åskådliggörs händelsen \(B=\{\min(X,Y)=2\}\). Notera att man använder den vektoriserade varianten av min, dvs pmin.

mutate(bivar1,B=pmin(X,Y)==2) %>%
  ggplot(aes(x=X,y=Y,fill=B,label=round(`p(x,y)`,2))) + 
    geom_raster(alpha=0.5) +
    geom_label()

Marginalfördelningar

Beräkning av marginalfördelningarnas massfunktioner \(p(x)\) och \(p(y)\) sker med hjälp av group_by() och summarise().

    bivar1px <-
      bivar1 %>%
      group_by(X) %>%
      summarise(`p(x)`=sum(`p(x,y)`))

    bivar1px %>%
      ggplot(aes(x=X,y=`p(x)`)) +
      geom_col(fill="blue") +
      labs(y="p(x)", title="Plot av marginalfördelningen p(x)")

Det här sättet att beräkna marginalfördelningen av \(X\) kan generaliseras till godtyckliga härledda variabler \(Z=h(X,Y)\). Här plottas fördelningen till variabeln \(Z=|X-Y|\).

    bivar1 %>% # Skapa fördelning av Z
      mutate(Z=abs(X-Y)) %>%
      group_by(Z) %>%
      summarise(`p(z)`=sum(`p(x,y)`)) %>% # plottning
      ggplot(aes(x=Z,y=`p(z)`)) +
      geom_col(fill="red") +
      labs(y="p(z)", title="Plot av marginalfördelning till Z=|X-Y|")

Betingade fördelningar

Man erhåller massfunktion \(p(y|x)=P(Y=y|X=x)\) för den diskreta betingade fördelningen \(P(Y|X)\) genom att dividera massfunktionen \(p(x,y)\) med marginalfördelningens \(p(x)\), dvs \[ p(y|x) = p(x,y)/p(x). \] Notera att man kan betrakta den betingade fördelningen \(P(Y|X)\) som en fördelning av \(Y\)-variabeln som beror på X-värdet. Man måste förstås kontrollera att \(p(x)\) inte är noll.

Vi utökar här tabellen bivar1 med hjälp av kommandot left_join(...) så att marginalfördelningen \(p(x)\) ingår som kolumn (med duplicerade värden). Sedan lägger vi lägger också till den betingade fördelningen \(p(y|x)\) som kolumn, med hjälp av mutate().

bivar1full <-
  bivar1 %>%
  left_join(bivar1px,by = "X") %>%
  group_by(X) %>%
  mutate(`p(x)`=sum(`p(x,y)`), `p(y|x)` = ifelse(`p(x)` > 0,`p(x,y)` / `p(x)`,0))
pr(bivar1full)

Betingade fördelningar

Här visas de betingade fördelningarna \(P(Y|X=x)\) för värdena \(x=0,1,2\) i en facet grid elelr facet_wrap().

ggplot(bivar1full,aes(x=Y,y=`p(y|x)`)) +
  geom_col(fill="blue") +
  facet_wrap(facets=vars(X)) +
  labs(title="De betingade fördelningarna av Y givet värden på X")

Uppgift 1

Med den diskreta fördelningen av \((X,Y)\) angiven ovan i bivar1.

ggplot(bivar2,aes(x=V,y=`p(v|u)`)) +
geom_col(fill="blue") +
geom_line()+
facet_wrap(facets=vars(U)) +
theme(panel.grid = element_line(color = "#8ccde3",
                                  size = 0.75,
                                  linetype = 2))+
labs(title="De betingade fördelningarna av V givet värden på U")
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

Uppgift 2

Kasta en tärning med värden i \(\{1, 2, . . . , 6\}\) ett första kast. Fortsätt sedan att kasta till du slår ett värde som är mindre eller lika med det första. Låt \(N\) ange antalet kast utöver det första och låt \(X\) ange värdet på första kastet och \(Y\) värdet på det sista kastet.

bivar12 <-  diceYN %>%
  group_by(Y) %>% 
  mutate(`p(E)`=(sum(N*(`p(n,y)`/`p(y)`) )))
sample_n(bivar12,4)
## # A tibble: 24 × 6
## # Groups:   Y [6]
##        Y     N    `p(n,y)` `p(y)`    `p(n|y)` `p(E)`
##    <int> <int>       <dbl>  <dbl>       <dbl>  <dbl>
##  1     1     3 0.0424       2.42  0.0175        3.41
##  2     1    13 0.00334      2.42  0.00138       3.41
##  3     1    12 0.00407      2.42  0.00168       3.41
##  4     1     8 0.00961      2.42  0.00396       3.41
##  5     2    14 0.000146     1.45  0.000101      2.03
##  6     2     3 0.0231       1.45  0.0160        2.03
##  7     2     3 0.0231       1.45  0.0160        2.03
##  8     2     9 0.00120      1.45  0.000825      2.03
##  9     3    16 0.000000850  0.950 0.000000894   1.52
## 10     3     2 0.0278       0.950 0.0292        1.52
## # … with 14 more rows
# Change dotsize and stack ratio
ggplot(bivar12, aes(x=`p(E)`, y=`p(n|y)`)) + 
  geom_dotplot(binaxis='y', stackdir='center',
               stackratio=1.8, dotsize=0.3, binwidth = 0.01)