Dette er en repetisjon av stoff fra innføringsemne i statistikk. Jeg har nedenfor lagt inn noen linjer R-kode, slik at du kan se hvordan vi kan bruke R til å gjøre noen enkle beregninger knyttet til binomiske sannsynligheter.

Tilfeldig variabel

En tilfeldig variabel er en kvantitativ størrelse som har et eller annet tilfeldig knyttet til sin verdi. Et av de første eksemplene man møter i statistikken er denne: La variabelen \(y\) være antall kron (eller mynt) som vi får ved å utføre \(n\) myntkast. Legg merke til at

Tilfeldige variable kan framkomme på mange andre måter også, men nå er det denne vi fokuserer på.

Binomisk variabel

Dersom alle myntkastene er uavhengige av hverandre og sannsynligheten for de to utfallene er identisk hver gang, så kalles \(y\) i en binomisk variabel (‘bi’ fordi det er 2 mulige utfall).

Hva mener vi med uavhengig? Jo, at utfallet av et myntkast ikke påvirker utfallet i de andre myntkastene. Tenk slik: Hvis jeg vet utfallet av mynkast 1, har jeg da noe informasjon om utfallet av myntkast 2? Hvis svaret er nei så er det uavhengighet.

Dersom vi gjør \(n=10\) myntkast, hva slags verdier kan \(y\) ha? Vel, vi kan risikere å få \(y=0\), selv om dette høres vanskelig ut (får bare mynt \(10\) ganger etter hverandre). Likeså, vi kan også få så stor verdi som \(y=10\), dersom alle myntkastene blir kron. Altså kan \(y\) få alle verdier \(0,1,2,...,10\). Men, ikke alle verdier er like sannsynlige! Dette beskrives av variabelens sannsynlighets-tetthet, og dersom vi plotter dette i et stolpediagram i R så ser det slik ut:

library(tidyverse)  # trenger ggplot2 pakken for å lage figur

binom.tbl <- tibble(antall_kron = 0:10) %>% 
  mutate(sannsynlighet = dbinom(x = antall_kron, size = 10, prob = 0.5))
ggplot(binom.tbl) +
  geom_col(aes(x = antall_kron, y = sannsynlighet)) +
  labs(x = "Antall kron", y = "Sannsynlighet") +
  scale_x_continuous(breaks = 0:10)

Til hver verdi som \(y\) kan ta er det tilordnet en sannsynlighet, og summen av alle sammen er \(1.0\). Formelen for denne tettheten finner du i STAT100-boka. Vi ser at \(y=5\) er det mest sannsynlige antallet, men at \(4\) og \(6\) er også vil opptre ganske ofte. Vi kan få fram tallene bak hver stolpehøyde ved å kikke i tabellen `binom.tbl:

flextable::regulartable(binom.tbl)

antall_kron

sannsynlighet

0

0.0009765625

1

0.0097656250

2

0.0439453125

3

0.1171875000

4

0.2050781250

5

0.2460937500

6

0.2050781250

7

0.1171875000

8

0.0439453125

9

0.0097656250

10

0.0009765625

Vi ser at sannsynligheten for at antall kron er 5 (tilsvarer \(y=5\)) har størst verdi, nesten \(0.25\). Vi kan også verifisere at disse tallene virkelig summerer seg til 1:

sum(binom.tbl$sannsynlighet)
## [1] 1

En binomisk variabel er altså antall utfall av den ene typen, der det er to typer (mulige utfall), og der vi gjentar dette \(n\) ganger. Sannsynligheten for den ene typen utfall trenger ikke være akkurat \(0.5\) (som for myntkast), men sannsynligheten for de to mulige utfallene må til sammen summere seg til \(1.0\).

Legg merke til at vi brukte funksjonen dbinom() for å få fram de binomiske sannsynlighetene. I R finner du funksjoner for å regne ut sannsynligheter for mange andre sannsynlighetsfordelinger også, som for eksempel normalfordelingen, eksponentialfordelingen, Poissonfordelingen osv. Sjekk ut ?dbinom for å lære om akkurat denne og liknende funksjoner som er nyttige i forbindelse med binomiske variable.

Forventningsverdi

En tilfeldig variabel vil (nesten) alltid ha en forventningsverdi. Tenk på dette som gjennomsnittsverdien vi får dersom vi har uendelig mange verdier for \(y\). Når vi har en sannsynlighets-tetthet som på figuren over, så tilsvarer forventningsverdien ‘tyngdepunktet’ i tettheten. Dersom du beregner et vektet gjennomsnitt av de \(11\) ulike verdiene som \(y\) kan ha, der hver verdi vektes med sin stolpehøyde, så får du eksakt forventningsverdien.

For den binomiske variabelen finnes det en enkel formel for å beregne forventningen:

\(E(y) = n\cdot p\)

der \(n\) er antall myntkast og \(p\) er sannsynligheten for kron hver gang. I figuren over er altså \(n=10\) og \(p=0.5\), som gir forventningsverdien \(5\).

Alt endrer seg med parameteren \(p\)

Sannsynligheten \(p\) kaller vi en parameter. Dette er sannsynligheten for den ene typen utfall. Den andre typen har da sannsynlighet \(1-p\). Denne bestemmer hvordan tettheten ser ut, og hva forventningsverdien blir. Her er et par binomiske tettheter med \(n=10\) som før, men der \(p\) er endret til \(0.25\) (til venstre) og \(0.90\) (til høyre).

binom1.tbl <- tibble(antall_kron = 0:10) %>% 
  mutate(sannsynlighet = dbinom(x = antall_kron, size = 10, prob = 0.25))
binom2.tbl <- tibble(antall_kron = 0:10) %>% 
  mutate(sannsynlighet = dbinom(x = antall_kron, size = 10, prob = 0.90))
bind_rows(binom1.tbl, binom2.tbl) %>% 
  mutate(p = c(rep("p=0.25", 11), rep("p=0.90", 11))) %>%
  ggplot() +
  geom_col(aes(x = antall_kron, y = sannsynlighet)) +
  labs(x = "Antall kron", y = "Sannsynlighet") +
  scale_x_continuous(breaks = 0:10) +
  facet_wrap(~p)

Legg merke til hvordan ‘toppen’ på fordelingene forskyver seg når sannsynligheten \(p\) blir mindre eller større enn \(0.5\). Forventningsverdien endrer seg også, og blir \(2.5\) og \(9.0\) for disse to tilfellene.

Sekvenseringsfeil

Vi har sekvensert reads som er \(300\) baser lange. Vi får oppgitt at feil-raten er 1%, dvs at sannsynligheten for feil på en gitt base er \(0.01\). For enhver base, så er det to muligheter: Enten er den feil eller ikke-feil. Vi antar også at i en read er basenes feil uavhengige, dvs. at dersom base 1 er feil, så ‘smitter’ ikke dette til base 2 eller 3 osv. Dette smaker av binomisk!

Vi ser på hver base som et myntkast, med sannsynlighet \(0.01\) for feil og \(0.99\) for ikke-feil. Antall feil på \(n=300\) baser er da en binomisk variabel \(y\). Forventningsverdien er da \(300\cdot 0.01 = 3\). Eller sagt på en annen måte: På en read av lengde \(n=300\) opptrer det i gjennomsnitt \(3\) sekvenseringsfeil.

Her er stolpediagrammet som viser tettheten for denne variabelen. Her har vi bare tatt med utfallene fra og med \(0\) til og med \(15\) ettersom utfallene \(16,17,...,300\) har så veldig små sannsynligheter.

sek_feil.tbl <- tibble(antall_feil = 0:15) %>% 
  mutate(sannsynlighet = dbinom(x = antall_feil, size = 300, prob = 0.01))
ggplot(sek_feil.tbl) +
  geom_col(aes(x = antall_feil, y = sannsynlighet)) +
  labs(x = "Antall sekvenseringsfeil i en read", y = "Sannsynlighet") +
  scale_x_continuous(breaks = 0:15)

Hva blir her sannsynligheten for at en read er feilfri? Dette er sannsynligheten for \(y=0\), som vi kan finne ved å se på tabellen sek_feil.tbl:

sek_feil.tbl %>% 
  filter(antall_feil == 0) %>% 
  flextable::regulartable()

antall_feil

sannsynlighet

0

0.04904089

Diskusjon

Dersom vi har Illumina-sekvenserte reads, så vil noen si at antall feil per read ikke beskrives godt av en binomisk variabel. Hva kan være grunnen til dette?

Vel, en forutsetning for at dette skal være binomisk er at sannsynligheten for feil er konstant langs hele read’en. Med Illumina vet vi at feilraten som regel er større mot slutten av hver read. Likevel gir binomisk variabel en brukbar tilnærming som vi anvender i stor stil.