c3.Rmd

---
title: "Rethinking. Поглавје 3. Белешки"
author: "discindo"
date: "January 20, 2020"
output:
  html_document:
    df_print: paged
    toc: yes
    toc_float: yes
    toc_depth: 2
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=TRUE)
```

```{r}
suppressPackageStartupMessages(library(rethinking))
suppressPackageStartupMessages(library(ggplot2))
```

# За

Целта на ова поглавје е да започнеме да размислуваме за моделирање во смисол на популации од параметри (примероци) наместо индивидуални тест статистики, и нормално, да научиме како да генерираме популации од проценки на параметри од постериорната дистрибуција. Главна причина за ова е дека моделите кои ни се потребни во вистинскиот свет се премногу комплицирани (многу параметри, со нерешени интеграли) за техниките како квадратна проценка или проценка со решетка. Оттука, потребно е да користиме симулатори (пр., Маркови вериги Монте Карло) кои во основа генерираат популации од параметри од постериорната дистрибуција. 

# Што направивме до сега?

Во претхнодото поглавје се запознавме со основите на веројатност и со Баесовото правило (теорема). Исто така разгледавме и една едноставна процедура за пресметка на постериорната веројатност кога моделот е едноставен. Примерот даден во книгата се однесува на биномен модел, како на пример фрлање паричка. Нѐ интересира која е веројатноста да добиеме 7 петки и 3 глави ако ја фрлиме фер паричката 10 пати? Фер паричка подразбира дека шансата да добиеме петка или глава е еднаква, 0.5 или 50%. Во `R`, и сите други програмски јазици, пресметување на оваа веројатност е едноставно. 

## Веројатност на петка или глава

```{r}
dbinom(
  x = 7,       # петки
  size = 10,   # вкупно фрлања
  prob = 0.5   # шансата за петка или глава на фер паричка
  )
```

ОК. Веројатноста да видиме 7 петки од 10 фрлања е `r dbinom(7, size=10, prob=0.5)`. Овие обзервации не се баш конзистентни со верна паричка. Ајде да провериме што ќе се случи доколку видиме 5 петки и 5 глави од 10 фрлања.

```{r}
dbinom(
  x = 5,       # петки  
  size = 10,   # вкупно фрлања
  prob = 0.5   # шансата за петка или глава на фер паричка
  )
```

Веројатноста во овој случај е поголема.

## Процена со решетка

Генерално, она што ќе не интересира е да ја најдеме најверојатната вредност на некој параметер откако сме собрале некои податоци. Ја фрламе паричката уште 90 пати и добиваме 67 петки и 33 глави Прашањето е дали е фер паричката? Или, погенерално, која е веројатноста на добивање петка имајќи ги во предвид видените податоци. Еден начин да го одговориме ова прашање е да тестираме неколку различни вредности за параметарот __веројатност на добивање петка__ и да видиме кој има најголема веројатност (кој е најможен).

```{r}
# вредности на параметарот кои сакаме да ги евалуираме
reshetka <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)

# повторување на пресметувањето за секој од можните вредности на параметарот
vozmoznost <- c()
for (i in seq_along(reshetka)) {
  vozmoznost[i] <- dbinom(67, 100, reshetka[i])
}
# сега правиме графикон од резултатите
plot(x = reshetka, y = vozmoznost, type = "b")
```

ОК. Веројатноста да ги видиме податоците 67 од 100 петки е најголема доколку веројатноста на добивање петка е околу 0.7 или 70%. Значи паричката најверојатно не е фер.

Горната процедура во која евауиравме неколку можни вредности за параметар што не интересира (најчесто еднакво распоредени во некој разумен распон) се нарекува grid approximation.

## Баесово ажурирање

Во Bayesian framework, оваа дистрибуција на веројатности треба да условиме со претходното верување (prior) и стандардизираме во однос на вкупната веројатнос на сите други настани (и за да добиеме веројатности помеѓу 0 и 1). Ајде да речеме дека немаме никава идеја за параметарот од интерес и да им дадеме иста претходна веројатност на сите вредности од решетката. Потоа, следејќи го баесовото правило, ја множиме веројатноста со претходните верувања и делиме со вкупната веројатност на сите можни настани. Стандардизацијата, покрај тоа што ни дава бројка помеѓу 0 и 1, ни кажува колкав дел од вкупната веројатност на сите настани отпаѓа на овој настан (параметар). Така може лесно да споредуваме вредности на параметри, некои ќе имаат постериорна веројатност од 10% други 1%.

```{r}
prethodni_verojatnosti <- rep(1, length(reshetka))
posteriorni_verojatnosti <- prethodni_verojatnosti * vozmoznost
posteriorni_verojatnosti <- posteriorni_verojatnosti / sum(posteriorni_verojatnosti)

par(mfrow = c(1,2))
plot(x = reshetka, y = prethodni_verojatnosti, type="b")
plot(x = reshetka, y = posteriorni_verojatnosti, type="b")
max(posteriorni_verojatnosti)
```

Ама ние имаме претходно искуство и знаеме дека некои вредности за претходната веројатност едноставно не се можни. На пример, добивме 67 петки и 33 глави, што значи веројатноста 0 за петка и 1 за глава (и обратно) не се можни. Имајќи ги овие информации, можеме да ја ажурираме низата на претходни веројатности (prior). Ајде да речеме дека вредности за веројатноста на глава под 0.2 и над 0.8 се невозможни:

```{r}
prethodni_verojatnosti <- c(0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0)
posteriorni_verojatnosti <- prethodni_verojatnosti * vozmoznost
posteriorni_verojatnosti <- posteriorni_verojatnosti / sum(posteriorni_verojatnosti)

par(mfrow=c(1,2))
plot(x=reshetka, y=prethodni_verojatnosti, type="b")
plot(x=reshetka, y=posteriorni_verojatnosti, type="b")

max(posteriorni_verojatnosti)
```

Забележете како податоците ја менуваат формата на кривата на веројатност. Тоа е Баесовото ажурирање (Bayesian updating).

Ако ја споредиме максималната постериорна веројатност со двата различни приори, ќе видиме дека приорот што беше поинформативен резултира со малку повисока постериорна веројатност. Ефектот тука е мал. За подобро да го осетиме овој ефект може да повториме пресметките со малку податоци.

```{r}
small_data <- c()
for (i in seq_along(reshetka)) {
  small_data[i] <- dbinom(6, 10, reshetka[i])
}
```

```{r}
prethodni_verojatnosti <- rep(1, length(reshetka))
posteriorni_verojatnosti <- prethodni_verojatnosti * small_data
posteriorni_verojatnosti <- posteriorni_verojatnosti / sum(posteriorni_verojatnosti)
plot(x=reshetka, y=posteriorni_verojatnosti, type="b")
max(posteriorni_verojatnosti)
```

```{r}
prethodni_verojatnosti <- c(0,0,0,1,1,1,1,1,0,0,0)
posteriorni_verojatnosti <- prethodni_verojatnosti * small_data
posteriorni_verojatnosti <- posteriorni_verojatnosti / sum(posteriorni_verojatnosti)
plot(x=reshetka, y=posteriorni_verojatnosti, type="b")
max(posteriorni_verojatnosti)
```

Промената во постериорната веројатност е позначителна и ефектот на претходното верување е евидентен во формата на кривата на постериори веројатност. 

Она што дополнително е интересно, е споредбата помеѓу разумен датасет (100 обзервации) и мал датасет (10 обзервации). Со големиот датасет, вредности за параметарот `< 0.5` имаат многу мала постериорна веројатност (0.0) иако претходната веројатност не беше нула! Додека со малку податоци, процедурата има помалку статистичка моќ, и дава постериорна веројатност дури и на вредности, како 0.3, кои интуитивно се поприлично неверојатни.

# Примероци од постериорната дистрибуција

Во примерот со глобусот, моделот е доволно едноставен така што можеме да ја проценуваме постериорната веројатност директно од проценката што ја добиваме од методот на решетка. Откако ја имаме идеја за формата на постериорната дистрибуција, земаме примерок од таа дистрибуција, пред сѐ, со цел да се запознаеме подетално со начинот на работа со вакви примероци, но и поради тоа што испитувањето на ваков начин ни дава одговори за постериорната дистрибуција кои што се важни за нашата работа / истражување.

Во иднина (т.е. со други методи), сѐ што ќе добиваме назад од проценката е некаков примерок од постериорната дистрибуција, но точните параметри на дистрибуцијата нема да ни бидат познати. Во многу случаи, зависно од нашиот приор и податоците, ваквите примероци од постериорната дистрибуција може да бидат прилично лоши или непрецизни. Дел од задачата тогаш ќе биде да утврдиме дека примерокот е соодветен и дека навистина се приближуваме до добри проценки на параметрите кои нѐ интересираат.

Кодот тука го следи моделот со глобусот од второто поглавје.

```{r}
# Креираме вектор можни вредности на параметарот за проценка со решетка
# Сакаме да ја најдеме вредноста која ја максимизира можноста
p_grid <- seq( from = 0 , to = 1 , length.out = 1000 )

# Користиме униформна приорна дистрибуција
# Вектор со илјада единици за секој елемент на мрежата
prior <- rep( x =  1 , times = 1000 )
```

Следно, ја пресметуваме возможноста на видените податоци (6 од 9 вода). За секоја вредност во решетката прашуваме, ако моделот е навистина биномен (претпоставка), тогаш колку е возможно да ги добиеме видените податоци, доколку пропорцијата на вода на земјата е некоја вредност (пр., 0.2).

```{r}
likelihood <- dbinom( x = 6 , size = 9 , prob = p_grid ) # 6 W; 3 L
```

Ја множиме возможноста со претходното верување (приор) за да го добиеме постериорот и стандардизираме со цел да имаме веројатност (вредност помеѓу 0 и 1). 

```{r}
posterior <- likelihood * prior 
# Стандардизираме
posterior <- posterior / sum(posterior) 
```

Постериорната веројатност ни кажува која е веројатноста на вредноста на параметарот (пропорција на вода на Земјата) имајќи во предвид видените податоци. 

Генерално, овој втор тип на веројатност е она што не интересира повеќе. Веројатноста на моделот при дадените податоци (Posterior probability) е поинтуитивен концепт отколку веројатноста на податоците при даден модел (Likelihood). 

График за постериорната дистрибуција пресметана според проценката на решетка.

```{r}
# plot( p_grid , 
#       posterior , 
#       type="b" ,
#       xlab="probability of water" , 
#       ylab="posterior probability" )

ggplot(data = data.frame(p_grid, posterior), aes(x = p_grid, y = posterior)) +
  geom_area(color = "steelblue", fill = "steelblue", alpha = .5)
```

Сега имаме постериорна дистрибуција. За овој едноставен модел, многу е лесно да земеме примерок бидејќи имаме предефинирани параметри (решетката) и нивната постериорна веројатност. Земаме примерок со замена, што значи истата вредност можеме да ја влкучиме повеќе пати, пропорционално со нејзината постериорна веројатност. Земаме 10.000 вредности од кои најчести се оние со најголема фрекфенција во постериорна дистрибуција.

```{r}
samples <- sample( p_grid,           # можните вредности на параметарот
                   prob = posterior, # веројантоста со која секоја вредност ќе биде извлечена
                   size = 1e4,       # колку вредности да извлечеме
                   replace = TRUE    # со замена 
                )

ggplot(data.frame(samples), aes(x = samples)) + 
  geom_histogram(color = "white", fill = "steelblue", alpha = 0.8, binwidth = 0.03)
# hist(samples) 
# dens(samples)
```


# Опишување на примероци

Штом ќе ја добиеме постериорната дистрибуцијата, моделот ја завршил работата. Но, сега треба да направи некакво толкување на дистрибуцијата, соодветно на потребите на истражувањето или работата поради која се правел моделот.

Генерално, има три групи на оценки на примероци, кои што даваат одговори на различен тип на прашања:
1. Интервали на граници;
2. Интервали на маса;
3. Проценки на точки.

## Граници

На пример, ако прашањето е колку е веројатноста пропорцијата на водата да е помала од 0.5, тоа лесно може да се пресмета од постериорната дистрибуција, која што за овој модел со глобус веќе ја имаме. 

```{r}
# add up posterior probability where p < 0.5
sum( posterior[ p_grid < 0.5 ] )
```

Со оглед на тоа дека дистрибуцијата вредности помеѓу 0 и 1 лесно е да пресметаме колку е веројатноста бидејќи во `R` функцијата `sum` ги собира вредностите на веројатноста каде што вредностите на параметарот `p` помал од `0.5`.

Но, како да ја добиеме веројатноста пропорцијата на водата да е помала од 0.5, кога би ја немале постериорната дистрибуција? Имаме 10.000 вредности од кои треба да видиме колку се помали од `0.5`. Правиме сума и делиме со 10.000 за да добиеме вредност помеѓу 0 и 1. 

```{r}
sum( samples < 0.5 ) / 1e4
```

## Маса

Интервали на маса се втората и, веројатно, варијантата со почеста употреба во науката. Барем таков е оној познатиот интервал што прикажува еднаква маса на веројатност на двете опашки и ако дистрибуцијата не е премногу асиметрична, дава добар приказ за формата на дистрибуцијата.

Така на пример ако сакаме да го видиме за под вредност на параметарот е дистрибуирана 80% од веројатноста на постериорната дистрибуција, можеме тоа да го добиеме со:

```{r}
quantile(samples , 0.8)
```

Пакетот `rethinking` што е придружен дел на овој курс содржи уште две готови функции кои што даваат интервали на маса и тоа PI (percetnile interval) и HDPI (highest density posterior interval). Првиот (PI) секогаш дава дво-опашест интервал така што во средината ја имаме масата на веројатноста што е еднаква на двете опашки заедно. 

На пример: 
```{r}
PI(samples, prob = 0.5)
```

Ни кажува дека 50% од веројатноста на параметарот (кој што во нашиот случај е пропорцијата на вода на Земјата) се наоѓа помеѓу `r quantile(samples, prob=0.25)` и `r quantile(samples, prob=0.75)`.

Додека HDPI ни кажува кој е намјмалиот интервал каде што се наоѓа наведената маса на дистрибуцијата. Во нашиот пример:

```{r}
(hpdi <- HPDI(samples, prob = 0.5))
```

Тоа е помеѓу `r hpdi[1]` и `r hpdi[2]`.

Генерално зборувано кога дистрибуцијата е симетрична нема голема разлика меѓу овие два интервали. Но, кога е лево или десно асиметрична (skewed) HDPI може да биде од поголема помош. За повеќе во врска со овие два интервали погледајте на стр. 57 од книгата.  

## Точки 

Последното резимирање на постериорната дистрибуција е доколку од неа се обидеме да извадиме точки, како што на пример се средната вредност, медијаната или модусот. Со оглед на тоа дека целиот пристап ни е да добиеме постериорна дистрибуција (а не точка) треба да имаме навистина добри причини зошто би сакале да се фокусираме на некоја точка.

Како и да е, ако веќе решиме, вредностите се добиваат лесно: 

```{r}
# mode of a bayesian sample (adj is a parameter passed to `density` internally)
rethinking::chainmode(samples, adj = 0.01) 
mean(samples)
median(samples)
```

# Адекватност на модели

Ако моделот што го користиме е соодветен за видените податоци, тогаш очекуваме дека доколку симулираме нови податоци врз база на нашиот модел и параметри, ќе добиеме симулирани податоци кои се многу слични на видените податоци. Овој факт може да го користиме за (1) да провериме дека кодот што сме го напишале го прави тоа што ние сакаме да го направиме и (2) да провериме дали моделот е адекватен за видените податоци.

Втората поента, адекватност на модел, е клучна. Честа пракса е да дефинираме неколку модели, да видиме које е нивната веројатност, и да ги рангираме од „најдобар“ до „најлош“. Оваа процедура е валидна во смисол на тоа дека ни кажува кој од моделите што сме ги земале во предвид е најдобар, или поточно, најмалку лош. Но не ни гарантира дека најмалку лошиот модел навистина ги опишува податоците добро, само дека ги опишува податоците подобро од другите модели што сме ги разгледале. За да провериме колку добро или лошо работи еден модел, можеме да симулираме податоци според тој модел и да ги споредиме со видените податоци. Доколку симулираните податоци се многу слични со видените, тогаш може да кажеме дека имаме адекватен модел. Постериорна предиктивна симулација (Posterior predictive simulation) е баесовата варијанта на оваа процедура. (Во frequentist likelihood framework оваа процедура се нарекува parametric bootstrap.)

Во пракса, оваа процедура е релативно едноставна. Земаме примерок од постериорната дистрибуција на параметрите на моделот (пропорционално со нивната постериорна веројатност, како и секогаш), и за секој сет на параметри во примерокот, симулираме одреден број на нови податоци. На пример може примарокот да е 1000 сетови на параметри, и за секој од нив да симулираме податоци еднаш или 1000 пати. Потоа ги опишуваме овие симулирани податоци заедно и сакаме да видиме колкав преклоп имаат со постериорната дистрибуција од видените податоци. Или пак калкулираме некоја карактеристика од симулираните сетови на податоци (популација од симулации) и проверуваме како, во просек, таа карактеристика личи споредбено со карактеристиката на видените податоци. Зависно од волуменот на податоци и комплексонста на моделот, оваа процедура може да бара многу временски и компјутерски ресурси, но имајќи ги во предвид клучните информации што ги добиваме, цената е прифатлива.

За биномниот модел е едноставно:

```{r}
# симулација за една можна вредност од постериорната дистрибуција
w <- rbinom(n = 1e4 , size = 9 , prob = 0.6 )
rethinking::simplehist(w)
```

Ако го погледнеме `w` ќе видиме бројки помеѓу 0 и 9, а не секвенци на единици и нули или W и L. Со оваа процедура симулираме _број на успеси од 9 обиди_, каде што веројатноста на успех (вода) е на пр. 0.6. Не ги симулираме примарите обзервации (фрлања на глобусот). 


```{r}
# симулација за вредностite во примерокот од постериорната дистрибуција
# (во овој пример немаме 1000 симулации за секоја вредност... )
w <- rbinom(n = 1e5 , size = 9, prob = samples)
simplehist(w)
```

За овој едноставен модел, процедурата е толку. За покомплицирани случаи, ќе треба да почекаме да завршат симулациите некое време.

Сега да ги опишеме симулираните податоци и да ги споредиме со видените податоци.


## Што ни кажуваат симулираните податоци

### Симулирање на број на успеси

```{r}
# geom_bar is more appropriate because the x axis is not continuous, values like 0.5 make no sense in a binomial model. You can't have 5.5 successes, has to be a whole number. If you use histogram, make sure the bin width is 1, so you don't get bars corresponding to fractions of 1.

ggplot(data.frame(w), aes(x=w)) + 
  # geom_histogram(binwidth = 1, fill="steelblue", color='white', alpha=.7) +
  geom_bar(fill="steelblue", width=.6, color="white", alpha=.7) +
  scale_x_continuous(breaks=0:9) +
  geom_vline(xintercept = 6, color="firebrick", size=5, alpha=.7)
```

Постериорниот предиктивен тест ни сугерира дека моделот врши добра работа. Дистрибуцијата на симулирани податоци (број на успеси) е центрирана на видената вредност (6 од 9), но со многу несигурност (можеме да видиме дека дистрибуцијата е раширена).

### Симулирање на примарни податоци

Но, сѐ што тестиравме е со која фрекфенција симулациите го вклучуваат бројот 6 (6 од 9 успеси). Во некои случаи ова е доволно. Ама честопати проблемот не е само фрекфенција на настан, туку, можеби секвенца на настани. 

На пример, во некоја интернет компанија прават тест на нов вебсајт, и гледаат дека стапката на кликање е 0.5, ама ги интересира дали секој втор корисник клика (секвенца 0101010101) или кликањата се групирани (може кликаат пријатели на фејсбук по препорака, и секвенцата е 1110011000). Прашањето е како да провериме дали симулираните податоци од биномниот модел ја следат шемата видените податоци. За овој тест ќе треба да симулираме примарни податоци (фрлања на глобусот), не бројот на успеси како погоре.

```{r}
voda <- 6/9

# за една вредност од постериорната дистрибуција
q <- sample(
  x = c("W", "L"),       # можните настани
  size = 9,              # број на фрлања на глобусот
  replace = TRUE,        # сакаме 9 елементи а имаме 2, мора со замена
  prob = c(voda, 1 - voda) # веројатност според видените податоци, 6 од 9 вода
  )

# за целата постериорна дистрибуција
z <- list()
for (i in seq_along(samples)) {
  voda <- samples[[i]]
  z[[i]] <- sample(
    x = c("W", "L"),       
    size = 9,              
    replace = TRUE,        
    prob = c(voda, 1 - voda) 
    )
}
```

Овој дел од кодот не од книгата. Во поглавјето не е даден код за овие симулации. Можеби има во дополнителен код во некои од repository-јата за курсот. Проверка е неопходна. Дали со овој код симулираме нешто слично на симулациите со `rbinom`?

```{r}
# прво броиме колку вода во секој симулиран датасет
count_water <- function(x) {sum(x == "W")}
count_water(z[[1]])
z_counts <- sapply(z, count_water)

water_counts <- data.frame(
  counts = c(w, z_counts), 
  type = c(rep("rbinom", length(w)), rep("sample", length(z_counts))))

ggplot(data = water_counts, aes(fill = type, x = counts)) + 
  geom_bar(color = "white", position = "dodge") +
  scale_x_continuous(breaks = 0:9)

# sanity check
# mean(w)
# sd(w)
# mean(z_counts)
# sd(z_counts)
# rethinking::HPDI(w)
# HPDI(z_counts)
```

ОК. Горе-долу добивме исти дистрибуциите на успеси, со тоа што имаме далеку помалки симулирани вредности (`z_counts`) отколку земени од постериорната дистрибуција (`w`). Значи со горниот код може да симулираме примарни обзервации на глобусот. Следи да го споредиме секвенците на "W" и "L" во симулираните податоци со видената секвенца.

```{r}
observed <- unlist(strsplit("WLWWWLWLW" , ''))
# ?rle

# Прво ја бараме должината на најдолгата секвенца од истиот елемент
# Значи за видените податоци очекуваме 3 (секвенцата WWW)
longest_stretch <- max(rle(observed)[[1]])

# Потоа бараме колку пати се сменил настанот
# Тоа е должината на резултатот од rle минус еден = 6
number_switches <- length(rle(observed)[[1]])-1
```

Следно, ја повторуваме оваа постапка за сите симулирани секвенци

```{r}
z_longest_stretch <- z_number_switches <- c()
for (i in seq_along(z)) {
  z_longest_stretch[[i]] <- max( rle( z[[i]] )[[1]] )
  z_number_switches[[i]] <- length( rle( z[[i]])[[1]] ) - 1
}
```

Ни останува само да ги опишеме овие дистрибуции и споредиме со видените вредности

```{r}
ggplot(data.frame(zls = z_longest_stretch), aes(x = zls)) + 
  geom_bar(width = 0.5, fill = "steelblue", alpha = 0.7) +
  geom_vline(xintercept = longest_stretch, size = 5, color = "firebrick", alpha = 0.7)
```

```{r}
ggplot(data.frame(zns=z_number_switches), aes(x=zns)) + 
  geom_bar(width = 0.5, fill="steelblue", alpha=.7) +
  geom_vline(xintercept = number_switches, size=5, color="firebrick", alpha=.7)
```

ОК. Слично како калкулациите во книгата, должината на повторување на настаните ја моделираме поприлично прецизно, но бројот на промени на настан во секвенцата не го моделираме адекватно.

Зависно од крајната цел на анализата, можеме да правиме препораки со овој модел, или да пробаме да го подобриме за да го моделира бројот на промени подобро.