03_inferential.Rmd

---
output:
  pdf_document: default
  html_document: default
---

# Järeldav statistika


```{r inferlibs}
library(tidyverse)
library(bayesboot)
library(boot)
library(rethinking) #HPDI() and PI(). Requires Stan.
#source("R/")
```
 
Kui EDA määrab graafiliste meetoditega andmete kvaliteeti ja püstitab uusi hüpoteese, siis järeldav statistika püüab formaalsete arvutuste abil vastata kahele lihtsale küsimusele: 1. mis võiks olla kõige usutavam parameetriväärtus? ja 2. kui suur ebakindlus seda hinnangut ümbritseb? Kuna andmed tulevad meile koos mõõtmisveaga ja koos bioloogilise varieeruvusega lõpliku suurusega valimina, on ebakindlus hinnagusse sisse ehitatud. Hea protseduur kvantifitseerib selle ebakindluse ausalt ja täpselt - siin ei ole eesmärk mitte niivõrd ebakindlust vähendada (seda tehakse eelkõige katse disaini tasemel), vaid seda kirjeldada.  
Järeldav statistika püüab, kasutades algoritme ja mudeleid, teha andmete põhjal järeldusi looduse kohta. 
Sellisel tegevusel on mõtet ainult siis, kui ühest küljest andmed peegeldavad tegelikkust ja teisest küljest tegelikkuses on rohkem, kui lihtsalt meie andmed. 
Kui andmed = tegelikkus, siis pole mõtet keerulisi mudeleid kasutada --- piisab lihtsast andmete kirjeldusest. 
Ja kui andmetel pole suurt midagi ühist tegelikkjsega, siis on need lihtsalt ebarelevatsed. 
Seega on järeldava statistika abil tehtud järeldused alati rohkem või vähem ebamäärased ning meil on vaja meetodit selle ebamäärasuse mõõtmiseks. Selle meetodi annab tõenäosusteooria.

## Järeldav statistika on tõenäosusteooria käepikendus

See õpik õpetab Bayesi statistikat, mis põhineb tõenäosusteoorial. Tänu sellele moodustab Bayesi statistika sidusa terviku, mille abil saab teha kõike seda, mida saab teha tõenäosusteooria abil. Bayesi statistika põhineb Bayesi teoreemil, mis on triviaalne tuletus tõenäosusteooria aksioomidest ja mille kohta on omakorda tõestatud teoreem, mis ütleb, et Bayesi teoreem on teoreetiliselt parim viis teha statistikat ja samal ajal mitte kaotada andmetes leiduvat informatsiooni. **NÕUAB TÄPSUSTUST!!!!!** 

Tõenäosusteooria on aksiomaatline süsteem, mille abil saame omistada numbriline väärtuse meie usu määrale mingisse hüpoteesi. Näiteks, kui me planeerime katset, kus me viskame kulli ja kirja ja teeme seda kaks korda, siis saame arvutada, millise tõenäosusega võime oodata katse tulemuseks kaht kirja. Aga seda tingimusel, et me võtame omaks mõned eeldused -- näiteks et münt on aus ja et need kaks viset on üksteisest sõltumatud. 

Sellel katsel on 4 võimalikku tulemust: H-H, H-T, T-H, T-T (H -kull, T - kiri). Tõenäosus saada 2-l mündiviskel 2 kirja, P(2 kirja) = 1/4, P(0 kirja) = 1/4 ja P(1 kiri) = 2/4 = 1/2. Sellega oleme andnud oma katseplaanile täieliku tõenäosusliku kirjelduse (pane tähele, et 1/4 + 1/4 + 1/2 = 1). Ükskõik kui keeruline on teie katseplaan, põhimõtteliselt käib selle analüüs samamoodi.
Tõenäosusteooria loomus seisneb kõikide võimalike sündmuste üleslugemises ning senikaua, kui me seda nüri järjekindlusega teeme, on vastus, mille me saame, tõsikindel.

Ehkki Bayesi statistika põhineb tõenäosusteoorial ja on sellega kooskõlas, ei ole see sama asi, mis tõenäosusteooria.
Statistikas pööratakse tõenäosusteoreetiline ülesanne pea peale ja küsitakse nii: kui me saime 2-l mündiviskel 2 kirja, siis millise tõenäosusega on münt aus (tasakaalus)? Erinevus tõenäosusteoreetilise ja statistilise lähenemise vahel seisneb selles, et kui tõenäosusteoorias me eeldame, et teame, kuidas süsteem on üles ehitatud, ja ennustame sellest lähtuvalt andmete tõenäosusi, siis statistikas me kontrollime neid eeldusi andmete põhjal. Seega annab tõenäosusteooria matemaatiliselt tõsikindlaid vastuseid ideaalmaailmade kohta, samas kui statistika püüab andmete põhjal teha kaudseid järeldusi päris maailma kohta. Selleks kasutame Bayesi teoreemi (vt allpool).

> Tõenäosusteooria määrab kõikide võimalike sündmuste esinemise tõenäosused, eeldades, et hüpotees H kehtib (H on siin lihtsalt teine nimi "eeldusele"). 

> Statistika arvutab H-i kehtimise tõenäosuse lähtuvalt kogutud andmetest, matemaatilistest mudelitest ning teaduslikest taustateadmistest.  

> Tõenäosus on defineeritud usu määrana mingisse hüpoteesi. Näiteks 62% tõenäosus (et populatsiooni keskmine pikkus < 180 cm) tähendab, et sa oled ratsionaalse olendina nõus kulutama mitte rohkem kui 62 senti kihlveo peale, mis võidu korral toob sulle sisse 1 EUR (seega 38 senti kasumit).  

**Tõenäosusteooria aksioomid** ütlevad tõlkes inimkeelde, et tõenäosused (P) jäävad 0 ja 1 vahele, et P(A) = 1 tähendab, et A on tõene, et P(A) = 0 tähendab, et A on väär, ning kui A ja B on ammendavad üksteist välistavad hüpoteesid, siis P(A) + P(B) = 1. Need aksioomid peaksid olema iseenesestmõistetavad ja ainult neist on tuletatud kogu tõenäosusteooria.

Need aksioomid, mis oma matemaatilises vormis postuleeriti Andrei Kolmogorovi poolt ca 1930, on tuletatavad järgmistest eeldustest:

+ ratsionaalne mõtlemine vastab kvalitatiivselt tervele mõistusele: lisatõendusmaterjal hüpoteesi kasuks tõstab selle hüpoteesi usutavust.

+ mõtlemine peab olema konsistentne: kui me võime järeldusi teha rohkem kui ühel viisil, peame lõpuks ikkagi alati samale lõppjäreldusele jõudma

+ kogu kättesaadav relevantne informatsioon tuleb järelduste tegemisel arvesse võtta (totaalse informatsiooni printsiip)

+ ekvivalentsed teadmised on representeeritud ekvivalentsete numbritega. 

Kui tõenäosused on 0 või 1, siis taandub tõenäosusteooria matemaatliselt oma erijuhule, milleks on lausearvutuslik loogika. Lausearvutuse oluline erinevus tõenäosusteooriast on, et kui selle abil on saavutatud valiidne tulemus, on see tõsikindel ja uued andmed ei saa põhimõtteliselt seda tulemust muuta. Seevastu tõenäosusteoorias ja statstikas muudavad uued andmed alati tõenäosusi. Selles mõttes ei saa tõenäosuslik teadus kunagi valmis.

### Tõenäosusteooria aksioomide formaalsed tulemid

1. Kui hüpotees A sisaldab endas hüpoteesi B, siis P(A) ei saa olla suurem kui P(B). 

2. kui A ja B on üksteist välistavad, siis tõenäosus, et A on tõene või B on tõene on nende kahe tõenäosuse summa --- *P(A v B) = P(A) + P(B)*. 

3. kui A ja B on üksteisest sõltumatud siis *P(A & B) = P(A) x P(B)*, kui A ja B ei ole üksteist välistavad, siis *P(A v B) = P(A) + P(B) – P(A & B)*. 

4. *P(A) = P(A & B) + P(A & mitte-B)* --- see on "totaalse tõenäosuse seadus" kitsendatuna kahele ammendavale teineteist välistavale hüpoteesile A & B, mille kohaselt me peame hüpotees A tõenäosuse leidmiseks arvestama ka B tõenäosusega.

5. *P(A | B) = P(A & B)/P(B)* kus *P(A | B)* on tinglik tõenäosus: A tõenäosus tingimusel, et B kehtib. *P(vihm | pilves ilm)* ei ole sama, mis *P(pilves ilm | vihm)*. 5. avaldisest saab lihtsa aritmeetika abil tuletada Bayesi teoreemi.

6. Bayesi teoreem: *P(A | B) = P(A)P(B | A)/P(B)*, kus *P(B) = P(A)P(B  A) + P(mitte-A)P(B  mitte-A)*  

Kui A on H~1~ ning mitte-A on ammendav ja välistav H~2~ ja B tähistab andmeid (data), saame Bayesi teoreemi ümber kirjutada 

*P(H~1~ | data) = ( P(H~1~)P(data | H~1~) )/( P(H~1~)P(data | H~1~) + P(H~2~)P(data | H~2~) )*

*P(H~1~ | data)* on *H~1~* kehtimise tõenäosus meie andmete korral, ehk posteerior, 

*P(H~1~)* on *H~1~* kehtimise eelnev, ehk meie andmetest sõltumatu, tõenäosus, ehk prior, 

*P(data | H~1~)* on andmete esinemise tõenäosus tingimusel, et H~1~ kehtib, ehk tõepära. 

Jagamistehe tehakse ainult selle pärast, et normaliseerida 1-le kõikide hüpoteeside tõenäosuste summa meie andmete korral ja seega viia posteerior vastavusse tõenäosusteooria aksioomidega --- kui meil on i ammendavat üksteist välistavat hüpoteesi, siis murrujoone alla läheb: *sum( P(data | H~i~)P(H~i~) ) = 1*. Bayesi teoreem on triviaalne tuletus tõenäosusteooria aksioomidest, milles pole mitte midagi maagilist. See ei ole automaatne meetod, mis tagaks inimkonna teadmiste kasvu, vaid lihtsalt parim võimalik viis andmemudeli ja taustateadmiste mudeli ühendamiseks ja normaliseerimiseks tinglikuks tõenäosuseks (hüpoteesi tõenäosus meie andmete ja taustateadmiste korral). Nüüd sõltub kõik mudelite, andmete ja taustateadmiste kvaliteedist. 

### Tõenäosusteooriast tulenevad statistika põhiprintsiibid

 1. statistilise analüüsi kvaliteet sõltub mudeli eeldustest & struktuurist. Kuna maailm ei koosne matemaatikast, teevad matemaatilised mudelid alati eeldusi maailma kohta, mis ei ole päris tõesed ja mida ei saa tingimata empiiriliselt kontrollida. Mündiviske näites eeldasime, et mündivisked olid üksteisest sõltumatud. Kui me sellest eeldusest loobume, läheb meie mudel keerulisemaks, sest me peame mudelisse lisama teavet visetevahelise korrelatsiooni kohta. Aga see keerulisem mudel toob sisse uued eeldused (vähemalt pool tosinat lisaeeldust). Üldiselt peaks mudeli struktuur kajastama katse struktuuri, mis kaasaegses statistikas tähendab sageli hierarhilisi mudeleid.
 
 2. statistilise analüüsi kvaliteet sõltub andmete hulgast. Kui kahe mündiviske asemel teeksime kakskümmend, siis saaksime samade eelduste põhjal teha oluliselt täpsemaid järeldusi mündi aususe kohta. 

 3. statistilise analüüsi kvaliteet sõltub andmete  kvaliteedist. Kui münt on aus, aga me viskame seda ebaausalt, siis, mida rohkem arv kordi me seda teeme, seda tugevamalt usub teadusüldsus selle tagajärjel millegisse, mis pole tõsi.

 4. statistilise analüüsi kvaliteet sõltub taustateadmiste kvaliteedist. Napid taustateadmised ei võimalda parandada andmete põhjal tehtud järeldusi juhul, kui andmed mingil põhjusel ei vasta tegelikkusele. Adekvaatsete taustateadmiste lisamine mudelisse aitab vältida mudelite üle-fittimist.

 5. Järeldused ühe hüpoteesi kohta mõjutavad järeldusi ka kõigi alternatiivsete hüpoteeside kohta. Relevantsete hüpoteeside eiramine viib ekslikele järeledustele kõigi teiste hüpoteeside kohta.  

## Andmed ei ole sama, mis tegelikkus


Kõigepealt illustreerime andmete ja tegelikkuse erinevust.
Simuleerimine on lahe sest simulatsioonid elavad mudeli väikeses maailmas, kus me teame 
täpselt, mida me teeme ja mida on selle tagajärjel oodata. Simulatsioonidega saame me 
hõlpsalt kontrollida, kas ja kuidas meie mudelid töötavad ning genereerida olukordi 
(parameetrite väärtuste kombinatsioone), mida suures maailmas kunagi ette ei tule. 
Selles mõttes on mudelid korraga nii väiksemad kui suuremad kui päris maailm. 

Alustuseks simuleerime juhuvalimi n=3 lõpmata suurest normaaljaotusega populatsioonist, mille keskmine on 100 ja sd on 20. Populatsioon (mis on statistiline mõiste) seisab siin tegelikkuse aseainena ja juhuvalim on meie andmete simulatsioon. Päris elus on korraliku juhuvalimi tõmbamine tehniliselt raske ettevõtmine ja, mis veelgi olulisem, me ei tea kunagi, milline on populatsiooni tõeline jaotus, keskmine ja sd. Elagu simulatsioon!

```{r}
set.seed(1) #makes random number generation reproducible
(Sample <- rnorm(n = 3, mean = 100, sd = 20)) #extra parentheses work as print()
mean(Sample)
sd(Sample)
```

Nagu näha on meie valimi keskmine 10% väiksem kui peaks ja valimi sd lausa kaks korda väiksem. Seega peegeldab meie valim halvasti populatsiooni --- aga me teame seda ainult tänu sellele, et tegu on simulatsiooniga.

Kui juba simuleerida, siis robinal: tõmbame ühe valimi asemel 10 000, arvutame seejärel 10 000 keskmist ja 10 000 sd-d ning vaatame nende statistikute jaotusi ja keskväärtusi. Simulatsioon on nagu tselluliit --- see on nii odav, et igaüks võib seda endale lubada. 

Meie lootus on, et kui meil on palju valimeid, millel kõigil on juhuslik viga, mis neid populatsiooni suhtes ühele või teisele poole kallutab, siis rohkem on valimeid, mis asuvad tõelisele populatsioonile pigem lähemal kui kaugemal. Samuti, kui valimiviga on juhuslik, siis satub umbkaudu sama palju valimeid tõelisest populatsiooniväärtusest ühele poole kui teisele poole ja vigade jaotus tuleb sümmeetriline. 

```{r}
N <- 3
N_simulations <- 10000
df <- tibble(a = rnorm(N * N_simulations, 100, 20), b = rep(1:N_simulations, each = N) )
Summary <-  df %>% group_by(b) %>% summarise(Mean=mean(a), SD= sd(a) ) 
Summary %>% ggplot(aes(Mean) ) + geom_histogram()
```

```{r}
mean(Summary$Mean) 
mean(Summary$SD)
```

Oh-hooo. Paljude valimite keskmiste keskmine ennustab väga täpselt populatsiooni keskmist aga sd-de keskmise keskmine alahindab populatsiooni sd-d. Valem, millega sd-d arvutatakse tõõtab lihtsalt kallutatult, kui n on väike (<10). Kui ei usu, korda eelnevat simulatsiooni valimiga, mille N=30.

Ja nüüd 10 000 SD keskväärtused:

```{r}
Summary %>% ggplot(aes(SD)) + geom_histogram()
```

```{r}
mode <-  function(x, adjust=1) {
  x <- na.omit(x)
  dx <- density(x, adjust=adjust)
y_max <- dx$x[which.max(dx$y)] 
y_max
}
mode(Summary$SD)
```

need koodiread seletavad lahti funtksiooni mode() koodi:
```{r}
dx <- density(rnorm(50)) #512 numbers x and y
head(dx$y)
plot(dx$x, dx$y)
plot(dx)
which.max(dx$y) #Determines the index of the maximum of a vector.
dx$x[which.max(dx$y)] 
```

SD-de jaotus on ebasümmeetriline ja mood ehk kõige tõenäolisem valimi sd väärtus, mida võiksime oodata, on u 14, samal ajal kui populatsiooni sd = 20. Lisaks on sd-de jaotusel paks saba, mis tagab, et tesest küljest pole ka vähetõenäoline, et meie valimi sd populatsiooni sd-d kõvasti üle hindab.

Arvutame, mitu % valimite sd-e keskmistest on > 25
```{r}
mean(Summary$SD>25)
```

Me saame >20% tõenäosusega pahasti ülehinnatud SD.

```{r}
mean(Summary$SD<15)
```

Ja me saame >40% tõenäosusega pahasti alahinnatud sd. Selline on väikeste valimite traagika.

Aga vähemalt populatsiooni keskmise saame me palju valimeid tõmmates ilusasti kätte --- ka väga väikeste valimitega.

Kahjuks pole meil ei vahendeid ega kannatust loodusest 10 000 valimi kogumiseks. Enamasti on meil üksainus valim. Õnneks pole sellest väga hullu, sest meil on olemas analoogne meetod, mis töötab üsna hästi ka ühe valimiga. Seda kutsutakse *bootstrappimiseks* ja selle võttis esimesena kasutusele parun von Münchausen. Too jutukas parun nimelt suutis end soomülkast iseenda patsi pidi välja tõmmata (koos hobusega), mis ongi bootstrappimise põhimõte. 
Statistika tõmbas oma saapaid pidi mülkast välja Brad Efron 1979 aastal.


```{r}
knitr::include_graphics("img/munchausen.jpg")
```


## Bootstrappimine

> Populatsioon on valimile sama, mis on valim bootstrappitud valimile. 

Nüüd alustame ühestainsast empiirilisest valimist ja genereerime sellest 1000 virtuaalset valimit. 
Selleks tõmbame me oma valimist virtuaalselt 1000 uut juhuvalimit (bootstrap valimit), millest igaüks on sama suur kui algne valim. 
Trikk seisneb selles, et bootstrap valimite tõmbamine käib asendusega, st iga empiirilise valimi element, mis bootstrap valimisse tõmmatakse, pannakse kohe algsesse valimisse tagasi. 
Seega saab seda elementi kohe uuesti samasse bootstrap valimisse tõmmata (kui juhus nii tahab). 
Seega sisaldab tüüpiline bootstrap valim osasid algse valimi numbreid mitmes korduses ja teisi üldse mitte. 
Iga bootstrap valimi põhjal arvutatakse meid huvitav statistik (näiteks keskväärtus) ja kõik need 1000 bootstrapitud statistikut plotitakse samamoodi, nagu me ennist tegime valimitega lõpmata suurest populatsioonist. 
Ainsad erinevused on, et bootstrapis võrdub andmekogu suurus, millest valimeid tõmmatakse, valimi suurusega ning, et iga bootstrapi valim on sama suur kui algne andmekogu (sest meie statistiku varieeruvus sõltub valimi suurusest ja me tahame seda varieeruvust oma bootstrapvalimiga tabada). 
Tüüpiliselt kasutatakse bootstrapitud statistikuid selleks, et arvutada usaldusintervall statistiku väärtusele.

> Bootstrap ei muuda meie hinnangut statistiku punktväärtusele. Ta annab hinnangu 
    ebakindluse määrale, mida me oma valimi põhjal peaksime selle punkthinnangu kohta tundma.
    
 Bootstrappimine on üldiselt väga hea meetod, mis sõltub väiksemast arvust eeldustest kui statistikas tavaks. Bootstrap ei eelda, et andmed on normaaljaotusega või mõne muu matemaatiliselt lihtsa jaotusega. Tema põhiline eeldus on, et valim peegeldab populatsiooni -- mis ei pruugi kehtida väikeste valimite korral ja kallutatud (mitte-juhuslike) valimite korral. Lisaks, tavaline bootstrap ei sobi hierarhiliste andmestruktuuride analüüsiks ega näiteks aegridade analüüsiks.

Bootstrap empiirilisele valimile suurusega n töötab nii:

1) tõmba empiirilisest valimist k uut virtuaalset valimit, igaüks suurusega n
2) arvuta keskmine, sd või mistahes muu statistik igale bootstrapi valimile. Tee seda k korda.
3) joonista oma statistiku väärtustest histogramm või density plot
4) nende andmete põhjal saab küsida palju toreidaid küsimusi --- vt allpool.

Mis on USA presidentide keskmine pikkus? Meil on valim 11 presidendi pikkusega.

```{r}
heights <- c(183, 192, 182, 183, 177, 185, 188, 188, 182, 185, 188) %>% as_tibble() 
n <- 11 #sample size
nr_boot_samples <- 1000 #the nr of bootstrap samples
a <- sample_n(heights, n * nr_boot_samples, replace=TRUE)  #create random sample with replacement
a$key <- rep(1:nr_boot_samples, each = n) #create a column "key" that cuts the sample into slices of size n
a1 <- a %>% group_by(key) %>% summarise(Value= mean(value)) #calculate the mean for each slice of n values
dens(a1$Value)
#HPDI(a1$Value)
```

Mida selline keskväärtuste jaotus tähendab? Me võime seda vaadelda posterioorse tõenäosusjaotusena. Selle tõlgenduse kohaselt iseloomustab see jaotus täpselt meie usku presidentide keskmise pikkuse kohta, niipalju kui see usk põhineb bootstrappimises kasutatud andmetel. Senikaua, kui meil pole muid relevantseid andmeid, on kõik, mida me usume teadvat USA presidentide keskmise pikkuse kohta, peidus selles jaotuses. Need pikkused, mille kohal jaotus on kõrgem, sisaldavad meie jaoks tõenäolisemalt tegelikku USA presidentide keskmist pikkust kui need pikkused, mille kohal posterioorne jaotus on madalam.   

Kuidas selle jaotusega edasi töötada? See on lihtne: meil on 1000 arvu ja me teeme nendega kõike seda, mida parasjagu tahame.

Näiteks me võime arvutada, millisesse pikkuste vahemikku jääb 92% meie usust USA presidentide tõelise keskmise pikkuse kohta. See tähendab, et teades seda vahemikku peaksime olema valmis maksma mitte rohkem kui 92 senti pileti eest, mis juhul kui USA presidentide keskmine pikkus tõesti jääb sinna vahemikku, toob meile võidu suuruses 1 EUR (ja 8 senti kasumit). Selline kihlveokontor on muide täiesti respektaabel ja akadeemiline tõenäosuse tõlgendus; see on paljude arvates lausa parim tõlgendus, mis meil on. 

Miks just 92% usaldusinterval? Vastus on, et miks mitte? Meil pole ühtegi universaalset põhjust eelistada üht usaldusvahemiku suurust teisele. Olgu meil usaldusinteval 90%, 92% või 95% --- tõlgendus on ikka sama. Nimelt, et me usume, et suure tõenäosusega jääb tegelik keskväärtus meie poolt arvutatud vahemikku. Mudeli ja maailma erinevused tingivad niikuinii selle, et konkreetne number ei kandu mudelist üle pärismaailma. NB! pane tähele, et eelnevalt mainitud kihlveokontor töötab mudeli maailmas, mitte teie kodu lähedasel hipodroomil.  

92% usaldusintervalli arvutamiseks on kaks meetodit, mis enamasti annavad vaid veidi erinevaid tulemusi.
1) HPDI --- Highest Density Probability Interval --- alustab jaotuse tipust (tippudest) ja isoleerib 92% jaotuse kõrgema(te) osa(de) pindalast
2) PI --- Probability Interval --- alustab jaotuse servadest ja isoleerib kummagist servast 4% jaotuse pindalast. See on sama, mis arvutada 4% ja 96% kvantiilid

```{r}
HPDI(a1$Value, prob=0.92) #highest probability density interval - non-symmetric
PI(a1$Value, prob=0.92) #symmetric method.
```


HPDI on üldiselt parem mõõdik kui PI, aga teatud juhtudel on seda raskem arvutada. Kui HPDI ja PI tugevalt erinevad, on hea mõte piirduda jaotuse enda avaldamisega --- jaotus ise sisaldab kogu informatsiooni, mis meil on oma statistiku väärtuse kohta. Intervallid on lihtsalt summaarsed statistikud andmete kokkuvõtlikuks esitamiseks.

Kui suure tõenäosusega on USA presidentide keskmine pikkus suurem kui USA populatsiooni meeste keskmine pikkus (178.3 cm mediaan)?
```{r}
mean(a1$Value>178.3)
```
Ligikaudu 100% tõenäosusega (valimis on 1 mees alla 182 cm, ja tema on 177 cm). Lühikesed jupatsid ei saa Ameerikamaal presidendiks!

### veidi keerulisem bootstrap

Eelnevalt tutvustasime nn percentiilmeetodit bootstrapi arvutamiseks. See lihtne meetod on küll populaarne ja annab sageli häid tulemusi, aga ei ole parim meetod bootstrappimiseks. Empiiriline bootstrap on sellest ainult veidike keerulisem aga see-eest annab robustsemaid tulemusi. Selles ei ploti me enam mitte 1000 statistiku väärtust vaid 1000 erinevust bootstrapitud statistiku väärtuse ja empiirilise valimi põhjal arvutatud statistiku väärtuse vahel.

```{r}
heights <- c(183, 192, 182, 183, 177, 185, 188, 188, 182, 185, 188) %>% as_tibble() 
n <- 11 #sample size
nr_boot_samples <- 1000 #the nr of bootstrap samples
a <- sample_n(heights, n * nr_boot_samples, replace=TRUE)  #create random sample with replacement
a$key <- rep(1:nr_boot_samples, each = n) #create a column "key" that cuts the sample into slices of size n
a1 <- a %>% group_by(key) %>% summarise(Value= mean(value) - mean(heights$value)) #calculate the difference between the mean for each slice of n values and the emprirical sample mean
dens(a1$Value)

```

Ja usaldusinterval tuleb niiviisi
```{r}
a <- HPDI(a1$Value, prob=0.95)
# a
# mean(heights$value)
CI <- c(mean(heights$value) - abs(a[1]), mean(heights$value) + abs(a[2]))
CI
```

### bayesboot()

See funktsioon pakub pisut moodsama meetodi --- Bayesian bootstrap --- mis töötab paremini väikeste valimite korral. Aga üldiselt on tulemused sarnased.
Hea lihtsa seletuse Bayesian bootstrapi kohta saab siit https://www.youtube.com/watch?v=WMAgzr99PKE ja lihtsa r koodi selle meetodi rakendamiseks saab siit https://www.r-bloggers.com/simple-bayesian-bootstrap/
Näited sellest, kuidas kasutada bayesbooti standardhälbe, korrelatsioonikoefitsiendi ja lineaarse mudeli koefitsientide usalduspiiride arvutamiseks leiate `?bayesboot` käsuga.

```{r}
library(bayesboot)
# Heights of the last 11 American presidents in cm (Kennedy to Trump).
heights <- c(183, 192, 182, 183, 177, 185, 188, 188, 182, 185, 188);
b1 <- bayesboot(heights, mean)
plot(b1, compVal=185)
#summary(b1)
#hist(b1)
HPDI(b1$V1, prob=0.95)
```


```{r, eval=FALSE}
# it's more efficient to use the a weighted statistic (but you can use a normal statistic like mean() or median() as well - as above).
b2 <- bayesboot(heights, weighted.mean, use.weights = TRUE)
```

Tõenäosus, et keskmine on suurem kui 182 cm
```{r}
# the probability that the mean is > 182 cm.
mean( b1[,1] > 182)
```

2 keskväärtuse erinevus (ES = keskmine1 - keskmine2):
 
```{r}
set.seed(1)
df <- tibble(a=rnorm(10), b=rnorm(10,1,1), c=a-b)
m <- bayesboot(df$c, weighted.mean, use.weights = TRUE )
plot(m, compVal=-1)

```


BayesianFirstAid raamatukogu funktsioon bayes.t.test() annab t-jaotuse mudeliga täpselt sama vastuse. See raamatukogu eeldab JAGS mcmc sämpleri installeerimist. Abi saab siit https://github.com/rasmusab/bayesian_first_aid ja siit https://faculty.washington.edu/jmiyamot/p548/installing.jags.pdf

```{r eval=FALSE}
library(BayesianFirstAid)
a <- bayes.t.test(df$a, df$b)
plot(a)
a1 <- a$mcmc_samples
a11 <- a1[[3]] %>% as.data.frame()
mean(a11$mu_diff > -1)
```

### Parameetriline bootstrap

Kui me arvame, et me teame, mis jaotusega on meie andmed, ja meil on suhteliselt vähe andmepunkte, võib olla mõistlik lisada bootstrapile andmete jaotuse mudel. Näiteks, meie USA presidentide pikkused võiksid olla umbkaudu normaaljaotusega (sest me teame, et USA meeste pikkused on seda). Seega fitime kõigepealt presidentide pikkusandmetega normaaljaotuse ja seejärel tõmbame bootsrap valimid sellest normaaljaotuse mudelist.
Normaaljaotuse mudelil on 2 parameetrit: keskmine (mu) ja standardhälve (sigma), mida saame fittida valimiandmete põhjal:
```{r}
mu <- mean(heights)
sigma <- sd(heights)
N <- length(heights)
sample_means <- tibble(value=rnorm(N*1000, mu, sigma), indeks=rep(1:1000, each=N)) 

sample_means1 <- sample_means %>% group_by(indeks) %>% summarise(Mean=mean(value))
    
ggplot(data = sample_means1, aes(x = Mean)) +
  geom_histogram(color = "white", bins = 20)
HPDI(sample_means1$Mean)
```

Üldiselt ei soovita me parameetrilist bootstrappi väga soojalt, sest täisbayesiaanlik alternatiiv, mida me kohe õppima asume, on sellest piisavalt paindlikum ja parem.

```{r, fig.cap="Bootstrap eesmärgiga leida usalduspiirid keskväärtusele",fig.align='center',out.width=600, echo=FALSE}
knitr::include_graphics("img/boot.pdf")
```

### Bootstrappimine ei ole kogu tõde

Bootstrappimine on võimas ja väga laia kasutusalaga meetodite kogum. Sellel on siiski üks oluline puudus. Nimelt arvestab bootstrap ainult andmetega ja ignoreerib taustateadmisi (parameetriline bootstrap küll eeldab taustateadmisetele tuginevalt jaotusmudelit, kuid ignoreerib kogu muud taustateadmist). Miks on see probleem?

Mõtleme hetkeks sellele teadusliku meetodi osale, millel põhineb suuresti näiteks Darwini liikide tekkimise argument. See on nn *inference to the best explanation*, mille kohaselt on eelistatud see teooria, mis on parimas kooskõlas faktidega, ehk mille kehtimise korral on meie andmete esinemine kõige tõenäolisem. Kui mõni hüpotees omistab andmete esinemisele suure tõenäosuse, siis me ütleme tehnilises keeles, et see hüpotees on **tõepärane** (*has high likelihood*). Esmapilgul tundub see kõik igati mõistlik, kuid proovime lihtsat mõtteeksperimenti. Selles juhtub nii, et loteriil võidab peaauhinna meile tundmatu kodanik Franz K. Meil on selle fakti (ehk nende andmete) seletamiseks kaks teooriat: 1. Franz K. võit oli juhuslik (loterii oli aus ja keegi peab ju võitma) ja 2. Franz K. vanem õde võltsis loterii tulemusi oma venna kasuks. Teine teooria sobib andmetega palju paremini kui esimene (sest kuigi keegi peab võitma, Franz K. võiduvõimalus oli väga väike); aga ometi eelistab enamus mõistlikke inimesi esimest teooriat. Põhjus on selles, et meil pole iseenesest mingit alust arvata, et Franz K.-l üldse on noorem õde, või et see õde omaks ligipääsu loteriile. Kui me aga saame teada, et Franz K. noorem õde tõesti korraldab loteriid, siis leiame kohe, et asi on kahtlane. 

Siit näeme, et lisaks tõepärale on selleks, et me usuksime mõne teooria kehtimisse, vaja veel, et see teooria oleks piisavalt tõenäoline meie taustateadmiste valguses. Bayesi teoreem ei tee muud, kui kvantifitseerib teooria kehtimise posterioorse tõenäosuse (järeltõenäosuse), kasutades selleks meie eelteadmiste ja tõepära kvantitatiivseid mudeleid. 
Seega, Bayesi paradigmas ei arvesta me mitte ainult andmetega, vaid ka taustateadmistega, sünteesides need kokku üheks posterioorseks jaotuseks ehk järeljaotuseks. Selle jaotuse arvutamine erineb bootstrapist, kuid tema tõlgendus ja praktiline töö sellega on samasugune. Erinevalt tavapärasest bootstrapist on Bayes parameetriline meetod, mis sõltub andmete modelleerimisest modeleerija poolt ette antud jaotustesse (normaaljaotus, t jaotus jne). Tegelikult peame me Bayesi arvutuseks modelleerima vähemalt kaks erinevat jaotust: andmete jaotus, mida me kutsume likelihoodiks ehk tõepäraks, ning eelneva teadmise mudel ehk prior, mida samuti modeleeritakse tõenäosusjaotusena.

> Bootstrapil (tavalisel ja Bayesi versioonil) on mõned imelikud formaalsed eeldused:
1. väärtused, mis ei esine valimiandmetes, on võimatud, 2. Väärtused, mis esinevad väljaspool valimi väärtuste vahemikku, on võimatud, 3. andmetes ei esine ajasõltuvusi ega hierarhilisi struktuure. Nendest puudustest hoolimata kasutatakse bootstrappimist laialt ja edukalt --- eelkõige tema lihtsuse ja paindlikuse tõttu. Küll aga tähendavad eelnimetatud puudused, et bootstrap on harva parim meetod teie ülesande lahendamiseks. 

Ehkki bootstrappimine ei arvesta taustateadmistega, ei tee seda olulisel määral ka paljud Bayesi mudelid (mudeldaja vaba valiku tõttu, mitte selle pärast, et mudel ei suudaks taustainfot inkorporeerida). Bayesi meetodite väljatöötajad ei tea sageli ette, milliste teaduslike probleemide lahendamiseks nende mudeleid hakatakse kasutama, ja seega ei kirjuta nad mudelisse ka väga ranget eelteadmist. Nende mudelite teadlastest kasutajad lepivad sageli selllega ja lasevad oma mudelite kaudu "andmetel kõneleda" enam-vähem sellistena, nagu need juhtuvad olema. Sellist lähenemist ei saa alati hukka mõista, sest vahest ei olegi meil palju eelteadmisi oma probleemi kohta, küll aga tuleb mainida, et sellistel juhtudel annab bootstrappimine sageli lihtsama vaevaga väga sarnase tulemuse, kui Bayesi täismäng. 

## Bayesi põhimõte

Bayesi arvutuseks on meil vaja teada

1) milline on "*parameter space*" ehk parameetriruum? Parameetriruum koosneb kõikidest loogiliselt võimalikest parameetriväärtustest. Näiteks kui me viskame ühe korra münti, koosneb parameetriruum kahest elemendist: 0 ja 1, ehk null kulli ja üks kull. See ammendab võimalike sündmuste nimekirja. Kui me aga hindame mõnd pidevat suurust (keskmine pikkus, tõenäosus 0 ja 1 vahel jms), koosneb parameetriruum lõpmata paljudest elementidest (arvudest). 

2) milline on "*likelihood function*" ehk tõepärafunktsioon? Me omistame igale parameetriruumi elemendile (igale võimalikule parameetri väärtusele) tõepära. Tõepära parameetri väärtusel x on tõenäosus, millega me võiksime kohata oma andmete keskväärtust, juhul kui x oleks see ainus päris õige parameetri väärtus. Teisisõnu, tõepära on kooskõla määr andmete ja parameetri väärtuse x vahel. Tõepära = Pr(andmed | parameetri väärtus). Näiteks, kui tõenäolised on meie andmed, kui USA keskmine president on juhtumisi 183.83629 cm pikkune? Kuna meil on vaja modeleerida tõepära igal võimalikul parameetri väärtusel (mida pideva suuruse puhul on lõpmatu hulk), siis kujutame tõepära  pideva funktsioonina (näiteks normaaljaotusena), mis täielikult katab parameetriruumi. Tõepärafunktsioon ei summeeru 100-le protsendile --- see on normaliseerimata.

3) milline on "*prior function*" ehk prior? Igale tõepära väärtusele peab vastama priori väärtus. Seega, kui tõepära on modelleeritud pideva funktsioonina, siis on ka prior pidev funktsioon (aga prior ei pea olema sama tüüpi funktsioon, kui tõepära). Erinevus tõepära ja priori vahel seisneb selles, et kui tõepärafunktsioon annab just meie andmete keskväärtuse tõenäosuse igal parameetriväärtusel, siis prior annab iga parameetriväärtuse tõenäosuse, sõltumata meie andmetest. See-eest arvestab prior kõikide teiste relevantsete andmetega, sünteesides taustateadmised ühte tõenäousmudelisse. Me omistame igale parameetriruumi väärtusele eelneva tõenäosuse, et see väärtus on üks ja ainus tõene väärtus. Prior jaotus summeerub 1-le. Prior kajastab meie konkreetsetest andmetest sõltumatut arvamust, kui suure tõenäosusega on just see parameetri väärtus tõene; seega seda, mida me usume enne oma andmete nägemist. Nendel parameetri väärtustel, kus prior (või tõepära) = 0%, on ka posteerior garanteeritult 0%. See tähendab, et kui te olete 100% kindel, et mingi sündmus on võimatu, siis ei suuda ka mäekõrgune hunnik uusi andmeid teie uskumust muuta (eelduselt, et te olete ratsionaalne inimene).

http://optics.eee.nottingham.ac.uk/match/uncertainty.php aitab praktikas priorit modelleerida (proovige *Roulette* meetodit).

Kui te eelnevast päriselt aru ei saanud, ärge muretsege. Varsti tulevad puust ja punaseks näited likelihoodi ja priori kohta.

Edasi on lihtne. Arvuti võtab tõepärafunktsiooni ja priori, korrutab need üksteisega läbi ning seejärel normaliseerib saadud jaotuse nii, et jaotusealune pindala võrdub ühega. Saadud tõenäosusjaotus ongi posteeriorne jaotus ehk posteerior ehk järeljaotus. Kogu lugu. 

Me teame juba pool sajandit, et Bayesi teoreem on sellisele ülesandele parim võimalik lahendus. Lihtsamad ülesanded lahendame me selle abil täiuslikult. Kuna parameetrite arvu kasvuga mudelis muutub Bayesi teoreemi läbiarvutamine eksponentsiaalselt arvutusmahukamaks (sest läbi tuleb arvutada mudeli kõikide parameetrite kõikide väärtuste kõikvõimalikud kombinatsioonid), oleme sunnitud vähegi keerulisemad ülesanded lahendama umbkaudu, asendades Bayesi teoreemi *ad hoc* MCMC algoritmiga, mis teie arvutis peituva propelleri Karlsoni kombel lendu saadab selle nimel, et tõmmata valim "otse" posterioorsest jaotusest. Meie poolt kasutatava MCMC *Hamiltonian Monte Carlo* mootori nimi on Stan. See on R-st eraldiseisev programm, millel on siiski R-i liides R-i pakettide rstan(), rethinking(), rstanarm() jt kaudu. Meie töötame ka edaspidi puhtalt R-s, mis automaatselt suunab meie mudelid ja muud andmed Stani, kus need läbi arvutatakse ja seejärel tulemused R-i tagasi saadetakse. Tulemuste töötlus ja graafiline esitus toimub jällegi R-is. Seega ei pea me ise kordagi Stani avama. 

Alustame siiski lihtsa näitega, mida saab käsitsi läbi arvutada.

### 1. näide

Me teame, et suremus haigusesse on 50% ja meil on palatis 3 patsienti, kes seda haigust põevad. Seega on meil kaks andmetükki (50% ja n=3). Küsimus: mitu meie patsienti oodatavalt hinge heidavad? Eeldusel, et meie patsiendid on iseseisvad (näiteks ei ole sugulased), on meil tüüpiline mündiviske olukord.

Parameetriruum: 0 surnud, 1 surnud, 2 surnud ja 3 surnud. Seega on meil neljaliikmeline parameetriruumi.

Edasi loeme üles kõik võimalikud sündmusteahelad, mis saavad loogiliselt juhtuda (et saada tõepärafunktsioon).

Me viskame kulli-kirja 3 korda: H - kiri, T - kull

Võimalikud sündmused on:

HHH,
HTH,
THH,
HHT,
HTT,
TTH,
THT,
TTT,



Kui Pr(H) = 0.5 ning H = elus ja T = surnud, siis lugedes kokku kõik võimalikud sündmused:

+ 0 surnud - 1,
+ 1 surnud - 3,
+ 2 surnud - 3,
+ 3 surnud - 1

Nüüd teame parameetriruumi igale liikme kohta, kui suure tõenäosusega me ootame selle realiseerumist. Näiteks, Pr(0 surnud) = 1/8, Pr(1 surnud) = 3/8, Pr(1 või 2 surnud) = 6/8 jne

Selle teadmise saab konverteerida tõepärafunktsiooniks

```{r}
x<-seq(from= 0, to=3) #parameter space as a grid
y<-c(1,3,3,1) #likelihood - 
plot(x,y, ylab="number of possibilities", xlab="number of deaths", type="b", main="likelihood")
```

Siit näeme, et üks surm ja kaks surma on sama tõenäolised ja üks surm on kolm korda tõenäolisem kui null surma (või kolm surma). 

Tõepära annab meile tõenäosuse Pr(mortality=0.5 & N=3) igale loogiliselt võimalikule surmade arvule (0 kuni 3). 

Me saame sama tulemuse kasutades binoomjaotuse mudelit (tegelikult kasutasime ka ennist sama mudelit). Ainus erinevus on, et nüüd on meil y teljel surmade tõenäosus.

```{r}
x<-seq(0, 3)
y <- dbinom(x,3, 0.5)
plot(x, y, type="b", xlab="nr of deaths", ylab="probability of x deaths", main="probability of x deaths out of 3 patients if Pr(Heads)=0.5")
```

Nüüd proovime seda koodi olukorras, kus meil on 9 patsienti ja suremus on 0.67

```{r}
x<-seq(0, 9)
y <- dbinom(x,9, 0.67)
plot(x, y, type="b", xlab="nr of deaths", ylab="probability of x deaths", main="probability of x out of 9 deaths if Pr(Heads)=0.67")
```

Lisame sellele tõepärafunktsioonile tasase priori (lihtsuse huvides) ja arvutame posterioorse jaotuse kasutades Bayesi teoreemi. Igale parameetri väärtusele on tõepära * prior proportsionaalne posterioorse tõenäosusega, et just see parameetri väärtus on see ainus tõene väärtus. Posterioorsed tõenäosused normaliseeritakse nii, et nad summeeruksid 1-le.

Me defineerime X telje kui rea 10-st arvust (0 kuni 9 surma) ja arvutame tõepära igale neist 10-st arvust. Sellega ammendame me kõik loogiliselt võimalikud parameetri väärtused.

```{r}
# define grid
p_grid <- seq( from=0 , to=9 )
# define flat prior
prior <- rep( 1 , 10 )
# compute likelihood at each value in grid
likelihood <- dbinom( p_grid , size=9 , prob=0.67 )

# compute product of likelihood and prior
unstd.posterior <- likelihood * prior 

# normalize the posterior, so that it sums to 1
posterior <- unstd.posterior/sum(unstd.posterior)
# sum(posterior) == 1 

plot( x = p_grid , y = posterior , type="b" ,
    xlab="nr of deaths" , ylab="posterior probability", main="posterior distribution" )

```

See on parim võimalik teadmine, mitu kirstu tasuks tellida, arvestades meie priori ja likelihoodi mudelitega. Näiteks, sedapalju, kui surmad ei ole üksteisest sõltumatud, on meie tõepäramudel (binoomjaotus) vale. 

### 2. näide: sõnastame oma probleemi ümber

Mis siis, kui me ei tea suremust ja tahaksime seda välja arvutada?
Kõik, mida me teame on, et 6 patsienti 9st surid.
Nüüd koosnevad andmed 9 patsiendi morbiidsusinfost (parameeter, mille väärtust me eelmises näites arvutasime) ja parameeter, mille väärtust me ei tea, on surmade üldine sagedus (see parameeter oli eelmises näites fikseeritud, ja seega kuulus andmete hulka).

Seega on meil
1) parameetriruum 0% kuni 100% suremus (0st 1-ni), mis sisaldab lõpmata palju numbreid.

2) kaks võimalikku sündmust (surnud, elus), seega binoomjaotusega modelleeritud tõepärafunktsioon. Nagu me juba teame, on r funktsioonis dbinom() kolm argumenti: surmade arv, patsientide koguarv ja surmade tõenäosus. Seekord oleme me fikseerinud esimesed kaks ja soovime arvutada kolmanda väärtuse.

3) tasane prior, mis ulatub 0 ja 1 vahel. Me valisime selle priori selleks, et mitte muuta tõepärafunktsiooni kuju. See ei tähenda, et me arvaksime, et tasane prior on mitteinformatiivne. Tasane prior tähendab, et me usume, et suremuse kõik väärtused 0 ja 1 vahel on võrdselt tõenäolised. See on vägagi informatsioonirohke (ebatavaline) viis maailma näha, ükskõik mis haiguse puhul!

**Tõepära parameetri väärtusel x on tõenäosus kohata meie andmeid, kui x on juhtumisi parameetri tegelik väärtus.**  Meie näites koosneb tõepärafunktsioon tõenäosustest, et kuus üheksast patsiendist surid igal võimalikul suremuse väärtusel (0...1). Kuna see on lõpmatu rida, teeme natuke sohki ja arvutame tõepära 20-l valitud suremuse väärtusel 

    Tehniliselt on sinu andmete tõepärafunktsioon agregeeritud iga üksiku 
    andmepunkti tõepärafunktsioonist. Seega vaatab Bayes igat andmepunkti 
    eraldi (andmete sisestamise järjekord ei loe).


```{r}
# define grid (mortality at 20 evenly spaced probabilities from 0 to 1)
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- rep( 1 , 20 )
# compute likelihood at each value in grid
likelihood <- dbinom( 6 , size=9 , prob=p_grid )
#plot prior
plot(p_grid, prior, type="b", main="prior")
#plot likelihood
plot(p_grid, likelihood, type="b", main="the likelihood function")
# compute product of likelihood and prior & standardize the posterior.
posterior <- likelihood * prior / sum(likelihood * prior)
#plot posterior
plot( x = p_grid , y = posterior , type="b" ,
    xlab="mortality" , ylab="posterior probability", main="posterior distribution" )
```

Nüüd on meil posterioorne tõenäosusfunktsioon, mis summeerub 1-le ja mis sisaldab kogu meie teadmist suremuse kohta.

> Alati on kasulik plottida kõik kolm funktsiooni (tõepära, prior ja posteerior).


#### Kui n=1

Bayes on lahe sest tema hinnangud väiksele N-le on loogiliselt sama pädevad kui suurele N-le. See ei ole nii klassikalises sageduslikus statistikas, kus paljud testid on välja töötatud N = Inf eeldusel ja töötavad halvasti väikeste valimitega.

Hea küll, me arvutame jälle suremust.

Bayes töötab andmepunkti kaupa (see et me talle ennist kõik andmed korraga ette andsime, on puhtalt meie mugavus)


```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- rep( 1 , 20 )
# compute likelihood at each value in grid
likelihood <- dbinom( 1 , size=1 , prob=p_grid )
posterior <- likelihood * prior / sum(likelihood * prior)

plot( x = p_grid , y = posterior , type="b" ,
    xlab="mortality" , ylab="posterior probability" )
```

esimene patsient suri - 0 mortaalsus ei ole enam loogiliselt võimalik (välja arvatud siis kui prior selle koha peal = 0) ja mortaalsus 100% on andmetega (tegelikult andmega) parimini kooskõlas. Posteerior on nulli ja 100% vahel sirge sest vähene sissepandud informatsioon lihtsalt ei võimalda enamat.


```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- posterior
# compute likelihood at each value in grid
likelihood <- dbinom( 1 , size=1 , prob=p_grid )
posterior1 <- likelihood * prior / sum(likelihood * prior)

plot( x = p_grid , y = posterior1 , type="b" ,
    xlab="mortality" , ylab="posterior probability" )
```

 Teine patsient suri. Nüüd ei ole 0 ja 1 vahel enam sirge posteerior. Posteerior on kaldu 100 protsendi poole, mis on ikka kõige tõenäolisem väärtus.


```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- posterior1
# compute likelihood at each value in grid
likelihood <- dbinom( 0 , size=1 , prob=p_grid )
# compute product of likelihood and prior
posterior2 <- likelihood * prior / sum(likelihood * prior)

plot( x = p_grid , y = posterior2 , type="b" ,
    xlab="mortality" , ylab="posterior probability" )
```

Kolmas patsient jäi ellu - 0 ja 100% mortaalsus on seega võimaluste nimekirjast maas ning suremus on ikka kaldu valimi keskmise poole (75%).

Teeme sedasama prioriga, mis ei ole tasane. See illustreerib tõsiasja, et kui N on väike siis domineerib prior posteerior jaotust. (Suure N korral on vastupidi, priori kuju on sageli vähetähtis.)


```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- c( seq(1:10), seq(from= 10, to= 1) )
# compute likelihood at each value in grid
likelihood <- dbinom( 1 , size=1 , prob=p_grid )
posterior <- likelihood * prior / sum(likelihood * prior)

plot(x=1:20, y=prior, type="b", main="prior")
plot( x = p_grid , y = posterior , type="b" ,
    xlab="mortality" , ylab="posterior probability", main="posterior")
```
1. patsient suri

```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- posterior
# compute likelihood at each value in grid
likelihood <- dbinom( 1 , size=1 , prob=p_grid )
# compute product of likelihood and prior
posterior1 <- likelihood * prior / sum(likelihood * prior)

plot(x=1:20, y=prior, type="b", main="prior")

plot( p_grid , posterior1 , type="b" ,
    xlab="mortality" , ylab="posterior probability", main="posterior" )
```

Teine patsient suri

```{r}
# define grid
p_grid <- seq( from=0 , to=1 , length.out=20 )
# define prior
prior <- posterior1
# compute likelihood at each value in grid
likelihood <- dbinom( 0 , size=1 , prob=p_grid )
posterior2 <- likelihood * prior / sum(likelihood * prior)

plot(x=1:20, y=prior, type="b", main="prior")
plot( x = p_grid , y = posterior2 , type="b" ,
    xlab="mortality" , ylab="posterior probability", main="posterior"  )
```
3. patsient jäi ellu. Nüüd on posteeriori tipp mitte 75% juures nagu ennist, vaid kuskil 50% juures --- tänu priorile.

## Mudelite keel

Siin vaatame kuidas kirjeldada mudelit nii, et masin selle ära tunneb. Meie mudelid töötavad läbi rethinking() paketi. See raamatukogu pakub kaks võimalust, kuidas mudelit arvutada, mis mõlemad kasutavad sama notatsiooni. Mõlemad võimalused arvutavad posteeriori mitte Bayesi teoreemi kasutades (nagu me ennist tegime), vaid kasutades stohhastilisi meetodeid, mis iseloomustavad posteeriori umbkaudu (aga piisavalt täpselt). Põhjuseks on, et keerulisemate mudelite korral on Bayesi teoreemi kasutamine liialt arvutusmahukas. 

Esiteks rethinking::map() leiab posteeriori tipu ja selle lähedal funktsiooni tõusunurga. Siin on eelduseks, et posteerior on normaaljaotus. See eeldus kehtib alati, kui nii prior kui tõepära on modelleeritud normaaljaotusena (ja ka paljudel muudel juhtudel). 

Teine võimalus on rethinking::map2stan(), mis suunab teie kirjutatud mudeli Stan-i. Stan teeb *Hamilonian Monte Carlo* simulatsiooni, kasutades valget maagiat selleks, et tõmmata valim otse posteerioorsest jaotusest. See on väga moodne lähenemine statistikale, töötab oluliselt aeglasemalt kui map, aga ei sõltu normaaljaotustest ning suudab arvutada hierarhilisi mudeleid, mis map-le üle jõu käivad. 

Me võime sama mudeli kirjelduse sõõta mõlemasse funktsiooni.

Lihtne mudel näeb välja niimodi:

   dead ~ dbinom(9, p) ,  # binomial likelihood
   
   p ~ dunif(0, 1)     # uniform prior
        

Tõepärafunktsioon on modeleeritud binoomjaotusena. Parameeter, mille väärtust määratakse on p, ehk suremus. See on ainus parameeter, mille väärtust me siin krutime.
NB! igale parameetrile peab vastama oma prior. Meil on selles mudelis täpselt 1 parameeter ja 1 prior. Vastuseks saame selle ainsa parameetri posterioorse jaotuse. Hiljem näeme, et kui meil on näiteks 452 parameetrit, mille väärtusi me koos arvutame, siis on meil ka 452 priorit ja 452 posterioorset jaotust.

```{r}
library(rethinking)
m1 <- map(
    alist(
        dead ~ dbinom(9, p) ,  # binomial likelihood
        p ~ dunif(0, 1)     # uniform prior
), data=list(dead=6) )

precis( m1 ) # summary of quadratic approximation
```


Nüüd tõmbame posteerioroorsest jaotusest valimi n=10 000. Selleks on funktsioon extract.samples()

```{r}
samples<-extract.samples(m1)
#hist(samples$p)
dens(samples$p)
#PI(samples$p, prob = 0.95) #leaves out equal 2.5% at both sides
HPDI(samples$p, prob = 0.95) #highest density 95% at the center
```

Kuus patsienti üheksast surid ja nüüd me usume, et tegelik suremus võib olla nii madal kui 37% ja nii kõrge kui 97%. Kui me tahame paremat hinnangut on meil vaja kas rohkem patsiente või informatiivsemat priorit (paremat taustainfot).


```{r}
library(rethinking)
m2 <- map(
    alist(
        dead ~ dbinom(90,p) ,  # binomial likelihood
        p ~ dunif(0,1)     # uniform prior
), data=list(dead=60) )
# display summary of quadratic approximation
precis( m2 )
```

```{r}
samples<-extract.samples(m2)
dens(samples$p)
#PI(samples$p, prob = 0.95) #leaves out equal 2.5% at both sides
HPDI(samples$p, prob = 0.95) #highest density 95% at the center
```

10 korda rohkem andmeid: nüüd on suremus määratud kuskile 57% ja 77% vahele (suure tõenäosusega)

## beta prior 

Nüüd anname sisse mõistlikuma struktuuriga priori: beta-jaotuse

Beta-prior katab vahemiku 0st 1ni ja sellel on 2 parameetrit, a ja b.

Siin mõned näited erinevatest beta parametriseeringutest
```{r}
x <- seq(0, 1, length = 1000)
plot(x, dbeta(x, 0.2, 0.2))
plot(x, dbeta(x, 1, 0.2))
plot(x, dbeta(x, 1, 1))
plot(x, dbeta(x, 2, 1))
plot(x, dbeta(x, 4, 1))
plot(x, dbeta(x, 2, 2))
plot(x, dbeta(x, 4, 4))
plot(x, dbeta(x, 200, 100))
```

 beta(θ | a, b) jaotuse keskväärtus on 
 
 *μ = a/(a + b)* 
 
ja mood on 

*ω = (a − 1)/(a + b − 2)* (juhul kui *a > 1* ja *b > 1*). 

Seega, kui a=b, siis on keskmine ja mood 0.5. Kui a > b, on keskmine ja mood > 0.5 ja kuid a < b, on mõlemad < 0.5. 

Beta jaotuse "laiuse" annab "konsentratsioon" *κ = a + b*. 
Mida suurem κ, seda kitsam jaotus.

*a = μκ*  

*b = (1 − μ)κ*

*a = ω(κ − 2) + 1* 

*b = (1 − ω)(κ − 2) + 1* for *κ > 2* 

Me võime κ-le omistada väärtuse nagu see oleks mündivisete arv, mis iseloomustab meie priori tugevust (juhul kui tõepära funktsioon tuleb andmetest, mis koosnevad selle sama mündi visetest). Kui meie jaoks piisaks ainult mõnest mündiviset, et priorist (eelnevast teadmisest) lahti ütelda, peaks meie prior sisaldama väikest kappat. 

    Näiteks, mu prior on, et münt on aus (μ = 0.5; a = b), aga ma ei ole selles väga kindel. 
    Niisiis ma arvan, et selle eelteadmise kaal võrdub sellega, kui ma oleksin näinud 8 mündiviske tulemust. 
    Seega κ = 8, mis tähendab, et a = μκ = 4 and b = (1 − μ)κ = 4. 
    Aga mis siis kui me tahame beta priorit, mille mood ω = 0.8 ja κ = 12? 
    Siis saame valemist, et a = 9 ja b = 3. 


```{r}
library(rethinking)
m3 <- map(
    alist(
        dead ~ dbinom(9,p) ,  # binomial likelihood
        p ~ dbeta(200,100)     # beta prior
), data=list(dead=6) )
# display summary of quadratic approximation
precis( m3 )

```


```{r}
samples<-extract.samples(m3)
dens(samples$p)
HPDI(samples$p, prob = 0.95) #highest density 95% at the center
```

Nagu näha on ka kitsa priori mõju üsna väika, isegi kui kui n=9.

## Lõpetuseks veel prioritest üldiselt. 

Neid võib jagada kolmeks: mitteinformatiivsed, väheinformatiivsed ehk "regularizing" ja informatiivsed. 

Mitteinformatiivseid prioreid ei ole sisuliselt olemas ja neid on soovitav vältida. Sageli kutsutakse tasaseid prioreid mitteinformatiivseteks. Neil on vähemalt 2 puudust. Tasane prior, mis ulatub lõpmatusse, on tehniliselt "improper", sest tema alune pindala ei summeeru ühele. Ja teiseks muudavad sellised priorid mcmc ahelad vähem efektiivseteks, mis võib teie arvutuse kihva keerata. 

Väheinformatiivsed priorid kujutavad endast kompromissi: nad muudavad võimalikult vähe tõepärafunktsiooni kuju, aga samas piiravad seda osa parameetriruumist, kust MCMC ahelad posteeriori otsivad (mis on soodne arvutuslikult). Nende priorite taga on filosoofiline eeldus, et teadlast huvitavad eelkõige tema enda andmed ja see, mida need ühe või teise hüpoteesi (parameetri väärtuse) kohta ütlevad. See eeldus on vaieldav, aga kui selle järgi käia, siis kulub vähem mõttejõudu eelteadmiste mudelisse formaliseerimiseks. 

Vähemalt suured farmaatsiafirmad seda hoiakut ei jaga ja kulutavad oma miljoneid informatiivsete priorite tootmiseks. Selles protsessis saavad kokku statistikud, teaduseksperdid ja psühholoogid, et inimkonna teadmisi võimalikult adekvaatselt vormida tõenäosusjaotustesse. Meie töötame siiski enamasti väheinformatiivsete prioritega.