05_rozdeleni.qmd

# Rozdělení pravděpodobnosti {#sec-rozdeleni}

::: callout-warning
## Cíle cvičení 

- Vědět kde v R najít rozdělení.
- Umět vysvětlit vztah mezi jednotlivými funkcemi rozdělení.
- Umět počítat s náhodnými čísly z rozdělení.
:::

V software R jsou všechna základní rozdělení součástí balíku `stats`. A jejich výčet získáme vyvoláním nápovědy `?distributions`. Pracovat s rozděleními budeme čtyřmi způsoby: - budeme generovat náhodná čísla z rozdělení - budeme pracovat s kvantilovou funkcí daného rozdělení - počítat s hustotou a pravděpododobnostní funkcí

Balík `stats` obsahuje mnoho rozdělení. **Alternativní**, **Binomické**, **Multinomické**, **Poissonovo**, **Normální**, **Fischer-Snedecorovo**, $\chi^2$, **Studentovo**, $t$, **Gumbellovo** a mnohá další. \
Podívejte se na `?Distributions`.

## Funkce rozdělení

Pro práci s rozděleními jsou k dispozici obecně tři funkce:

```{r eval=FALSE}
r___() # generování náhodných čísel z rozdělení
d___() # funkce hustoty rozdělení
p___() # Pravděpodobnostní funkce 
q___() # Kvantilová funkce rozdělení
```


Pro úplnost jsou zde uvedeny funkce pro permutace `factorial()` a kombinace `choose()`.

```{r,fig.align='center'}
par(mfrow = c(1,2))
plot(x = 0:6, 
     y = factorial(0:6), # <1>
     type = "s",
     main = "Permutace", 
     ylab = bquote("log(0!:6!)"), 
     xlab = "",
     log = "y")
plot(x = 0:6, 
     y = choose(6, 0:6), # <2>
     type = "s", 
     xlab = "",
     main = "Binomické koeficienty")
```
1. Kolika způsoby je možné seřadit $n\leq6$ prvků?
2. Kolika způsoby lze vybrat $n-$tici z $6$ prvků?

### Distribuční funkce rozdělení

Rozdělení pravděpodobnosti náhodné veličiny lze jednoznačně popsat tzv. distribuční funkcí @eq-distribucni-funkce.

$$
F(x) = P(X\leq x) = P(\omega_i\in\Omega : X(\omega_i)\leq x)
$$ {#eq-distribucni-funkce}

#### Empirická distribuční funkce {#sec-rozdeleni-empiricka}

Sestavuje se na základě sesbíraného vzorku dat, sděluje jak velký podíl dat je stejných
nebo menších než konkrétní hodnota. K jejímu vytvoření potřebujeme data seřazená.

```{r, fig.align='center', fig.width=10}
Qd <- read.fwf(file = "./data/01138000.dly", #<1>
               widths = c(8, rep(10, times = 5)))[, 5] #<1>
```

$$
P(X\leq x_k) \approx \dfrac{k - 0.3}{n + 0.4}
$$
kde $k$ je pořadí v seřazeném souboru a $n$ je počet prvků v souboru.

::: callout-tip
## Úloha

1.  a) Nahrajte data. \     
    b) Vytvořte vektor pořadí hodnot. \  
    c) Spočítejte pravděpodobnosti. \  
    d) Vyneste do grafu. \   
    d) Proveďte s pomocí funkce `ecdf()`.
:::


```{r, fig.align='center', fig.width=10, echo=FALSE}
k <- rank(Qd)
n <- length(Qd)
P <- (k - 0.3)/(n + 0.4)

dat <- data.frame(k, n, P, Qd)
dat <- dat[order(dat$P), ]
par(mfrow = c(1, 2))
plot(y = dat$P, # <2>
     x = dat$Qd, # <2>
     ylab = expression(P(X<=x[k])), # <2>
     xlab = expression(T[avg]),  # <2>
     type = "l",
     lwd = 1.25) # <2>

plot(x = dat$Qd, # <2>
     y = 1-dat$P, # <2>
     ylab = expression(P(X>x[k])), # <2>
     xlab = expression(T[avg]),  # <2>
     type = "l", 
     lwd = 1.25) # <2>

```



### Kvantilová funkce

Kvantilová funkce je inverzní funkcí k distribuční. Pokud distribuční funkce
udává s jakou pravděpodobností bude hodnota náhodného pokusu menší nebo rovna $x$,
kvantilová funkce udává, pro jaké $x$ bude výsledek náhodného pokusu s danou pravděpodobností
$y$ měnší nebo roven $x$.

```{r}
qnorm(p = 0.8, mean = 2, sd = 1) #<1>
```
1. S $80\%$ pravděpodobností bude hodnota $\leq$`r qnorm(p = 0.8, mean = 2, sd = 1)`.

### Náhodná čísla z rozdělení

Pro generování náhodných čísel lze použít rozdělení. 

```{r,eval=FALSE}
runif(n = 10, min = 0, max = 1) # <1>
rpois(n = 15, lambda = 2.4) # <2>
```
1. Generování 10 čísel z rovnoměrného rozdělení z intervalu $(0;1)$
2. Generování 15 čísel z Poisonova rozdělení z intervalu $(0;1)$

Generovaná čísla nejsou náhodná v pravém slova smyslu, ale označují se jako *pseudonáhodná*, neboť při jejich tvorbě se vychází z jiste sekvence čísel. Tuto sekvenci je možné přímo zvolit, čímž je zajištěna 
funkcí `set.seed()` volanou před každou vygenerovanou sekvencí.

```{r}
x <- sample(x = 1:1e3, size = 1) # <1>
set.seed(x)  # <1>
runif(1)

set.seed(x)  # <2>
runif(2)
```


```{r}
seed <- .Random.seed
head(seed, 10)
```
```{r,fig.align='center'}
plot(.Random.seed)
```


## Diskrétní rozdělení

<!-- ### Alternativní rozdělení -->
<!-- ### Rovnoměrné rozdělení -->
<!-- ### Binomické rozdělení -->
### Poissonovo rozdělení

```{r, fig.align='center'}
# par(mfrow = c(2, 2))
lambda <- c(1, 4, 8)
plot(x = 1:20, ylim = c(0,0.4), type = "n")
for(l in lambda) {
  lines(dpois(x = 1:20, lambda = l), lty = l)
  points(dpois(x = 1:20, lambda = l), col = l)
}
```


::: callout-tip
## Úloha

  8. Vytvořte jednoduchý simulátor deště pro jeden rok podle následujícího zadání\
    a) Buď prší, nebo ne. Prvděpodobnost že neprší, je 90 \%.\
    b) Pokud pršelo předchozí časový úsek, tak se zvyšuje pravděpodobnost že prší 
  dnes na 25 \%.
:::

Řešení:\
```{r, fig.align='center'}
#| code-fold: true
x <- vector(mode = "numeric", length = 365) #<1>

for(i in seq_along(x)) {
  if(i == 1) next #<2> 
  x[i] <- ifelse(
            test = x[i - 1] > 0, #<3>
            yes = rbinom(n = 1, #<3>
                         size = 1, #<3>
                         prob = 0.10) * rchisq(n = 1, #<3>
                                               df = 1, #<3>
                                               ncp = 0), #<3>
            no = rbinom(n = 1, 
                        size = 1, 
                        prob = 0.25) * rchisq(n = 1, 
                                              df = 3, 
                                              ncp = 1)) 
}
barplot(x, 
        ylim = 1.5*c(max(x), 0), 
        col = "dodgerblue4", 
        border = "dodgerblue4")
```
1. Předem alokuejeme vektor hodnot.
2. Je třeba ošetřit případ, kdy začínáme a neexistuje úsek $x(t-1)$.
3. Pokud pršelo v předchozím čase



## Spojitá rozdělení

Během počítání s rozděleními se snažíme nalézt odpovědi na podobné otázky:\
"Jaká je pravděpodobnost, že veličina $X$ je menší než 7.5?"

### Normální rozdělení

Chceme spočítat pravědpodobnost, že náhodná hodnota, vybraná z populace $X\sim \mathsf{N}(\mu = 3; \sigma^2 = 2,2)$ bude menší než 5,5. Početně je úkon řešen standardizací a určením nového 

::: {layout-ncol=2}
```{r,echo=FALSE, fig.cap="Pravděpodobnost, že $x$ je menší nebo rovno 2."}
x_values <- seq(from = -5+3, to = 2+3, by = 0.001)
DNORM <- dnorm(x_values)

curve(dnorm(x, 3, sqrt(2.2)), -5+3, 2+3, 
      main = "Normální rozdělení", 
      cex.main = 1.5, #expression(paste(mu == 3, "; ", sigma^2 == 2,2))
      cex.axis = 1.5,
      adj = 0, 
      xlab = "", 
      ylab = "")

polygon(x = c(x_values, 2+3, -5+3),
        y = c(DNORM, 3, sqrt(2.2)),    
        col = "gray95", 
        border = 1)
segments(x0 = 2, x1 = 2, y0 = -0.2, y1 = dnorm(5.5))
arrows(x0 = 2, x1 = -2.2, y0 = 0.02, y1 = 0.02, length = 0.06, lwd = 1.2, code = 2, lty = 1)
text(x = 2.5, y = 0.25, expression(paste(X%~%N("3; 2,2"))), cex = 1.5)
```

```{r, echo=FALSE, fig.cap="Pravděpodobnost, že $x$ je menší nebo rovno 2."}
x_values <- seq(from = -5, to = 2, by = 0.001)
DNORM <- dnorm(x_values)

curve(dnorm(x), -5, 5, 
      main = "Standardizované Normální rozdělení", 
      cex.main = 1.5, #expression(paste(mu == 0, "; ", sigma == 1)) 
      cex.axis = 1.5,
      adj = 0, 
      xlab = "", 
      ylab = "")

polygon(x = c(x_values, 2, -5),
        y = c(DNORM, 0, 0),    
        col = "gray95", 
        border = 1)
segments(x0 = 2, x1 = 2, y0 = -0.2, y1 = dnorm(2))
arrows(x0 = 2, x1 = -2.2, y0 = 0.02, y1 = 0.02, length = 0.06, lwd = 1.2, code = 2, lty = 1)
text(x = 2.5, y = 0.25, expression(paste(X%~%N("0; 1"))), cex = 1.5)
```
:::

Provedeme standardizaci na $Z\sim\mathsf{N}(0; 1)$
$$
Z = \dfrac{X - \mu}{\sigma} \approx \dfrac{5,5 - 3}{1,48324} \approx 1,6855
$$
$P(Z < 1,69)$ Spočítáme s pomocí kvantilu Standardizovaného Normálního rozdělení.
Tento úkon R dělá z nás, ale pro porovnání si uvedeme obě varianty:
```{r}
cbind(
  N = pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE),
  Z = pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE)
)
```

Hledaná pravděpodobnost je `r round(pnorm(q = 1.69, mean = 0, sd = 1, lower.tail = TRUE)*100, 2)` \%.


::: callout-note
## Poznámka

Zde je dobré upozornit na numerické zaokrouhlovací chyby, které během počítačového zpracování dat 
vznikají. 
```{r}
all.equal(
  pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE), 
  pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE)
)
all.equal(
  pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE), 
  pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE), tolerance = 0.005
)
```
:::





::: callout-tip
## Úloha

1.  a) S pomocí funkce `curve()` vykreslete hustotu funkce normovaného normálního rozdělení
danou předpisem:
$$
f(x) = \dfrac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}
$$
    b) Máme náhodnou veličinu $X$, jež sleduje Normální rozdělení s parametry $\mu = 5$ a $\sigma^2 = 4,13$, nebo-li $X\sim \mathsf{N}(5; 4,13)$. Určete pravděpodobnost:\     
        - $P(X < 8)=$ `r fitb(c("0.7662008", "0,7662008"))`\         
        - $P(X \leq 8)=$ `r fitb(c("0.7662008", "0,7662008"))`\         
        - $P(X > 12,4)=$ `r fitb(c("0.0365851", "0,0365851"))`\      
    c) Pro jakou hodnotu $\gamma$ je pravděpodobnost přibližně:
        - $P(X < \gamma) = 0,76\quad\gamma=$ `r fitb(c("7.92", "7,92"))`\          
        - $P(X > \gamma) = 0,16\quad\gamma=$ `r fitb(c("9.10", "9,10"))`            
:::

<!-- Řešení: -->
<!-- ```{r} -->
<!-- #| code-fold: true -->
<!-- curve(expr = 1/(sqrt(2*pi))*exp(-(x^2)/(2*1^2)), from = -3, to = 3) -->
<!-- curve(expr = 1/(sqrt(2*pi))*exp(-(x^2)/(2*1^2)), from = -3, to = 3) -->
<!-- ``` -->

### Studentovo-$t$ rozdělení

Používá se pro odhad střední hodnoty souboru pocházejícího z normálního rozdělení 
s neznámým parametrem $\sigma$. Více v @sec-testovani.

```{r}
qt(p = 0.95, df = 100)
```


### Logaritmicko-normální rozdělení

Nabývá pouze kladných hodnot. Můžeme jej nalézt například v rozdělení **průtoků, 
simulacích chemické konncentrace.


::: callout-note
## Poznámka

ČHMÚ používá pro odhad a konstrukci $m-$denních průtoků logaritmicko-normální 
rodzdělení s pěti parametry.     \

tvar transformace se třemi, s pěti a ve standardizovaném tvaru pro průtoky
<!-- $$ -->
<!-- y = e^x + y_o\\ -->
<!-- y = a\cdot e^{sign(x)\cdot|x|^b}+y_0\\ -->
<!-- q_{norm} = \alpha\cdot e^{sign(x)\cdot|x|^\beta}+q_{100\%} -->
<!-- $$ -->

  - $\mu$: výběrový průmer normálního rozdělení   \
  - $\sigma$: směrodatná odchylka normálního rozdělení   \
  - $\gamma$: parametru posunu   \
  - $\alpha$: parametr měřítka   \
  - $\beta$: parametr tvaru (3. a 4. moment)\

$$
f(x;\mu, \sigma, \gamma, \alpha, \beta) = \dfrac{\beta}{(x-\gamma)\sigma\sqrt{2\pi}}\exp\left(-\frac{\left(\ln{\left(\dfrac{x-\gamma}{\alpha}\right)}-\mu\right)^2}{2\sigma^2}\right),\quad x>\gamma
$$
```{r, echo=FALSE, fig.align='center', }
# x <- 1:100
ln5 <- function(x) ((beta / ((x - gamma) * sigma * sqrt(2 * pi))) * exp(-((log((x - gamma) / alpha) - mu)^2) / (2 * sigma^2)))

mu <- 0
sigma <- 1
gamma <- 0.8
alpha <- 1
beta <- 2
curve(expr = ln5, 
      from = gamma + 0.01, 
      to = 20, 
      n = 1000, 
      xname = "x", 
      col = "#102d3d", 
      main = "Hustota funkce logaritmicko-normálního rozdělení s 5 parametry", 
      xlab = "", 
      ylab = "f(x)")
```


:::


### $F$ rozdělení
Jinak také *Fisherovo-Snedecorovo* rozdělení je používáno pro sestrojení $100(1-\alpha)\%$ 
intevalu spolehlivosti pro podíl rozptylů normálního rozdělení a je to modelové rozdělení
testovací statistiky pro ověření shodnosti dvou rozptylů. V @sec-anova nalezneme $F-$rozdělení 
u testování hypotézy o rovnosti středních hodnot u více než dvou výběrových souborů.

### $\chi^2$ rozdělení

Narozdíl od předchozích zde uvedených rozdělení není $\chi^2$ rozdělení u pozorovaných veličin příliš časté. 
Nicméně je to významné rozdělení z hlediska testování statistických hypotéz (@sec-testovani) při porovnávání
rozdělení vzájemně. Používá se při stanovení intervalů spolehlivosti výběrový rozptyl.

$$
\dfrac{(n-1)s_x^2}{\chi^2_{1-\alpha/2}} \leq \sigma_x^2\leq \dfrac{(n-1)s_x^2}{\chi^2_{\alpha/2}}
$$ {#eq-odhad-chisq-rozptyl}

::: callout-tip
## Cvičení
  1.  Generujte 10 čísel s pomocí normálního rozdělení s parametry $\mu = 2.3$ a $\sigma = 4.23$.\
  2.  Generujte 50 čísel $X\sim\mathsf{Po}(\lambda = 2.3)$ s pomocí `set.seed(123)`. Doplňte střední hodnotu 
      generovaného souboru `r fitb(2.36)`    \    
  3.  Jaká je pravděpodobnost, že veličina $X\sim\mathsf{N}(1.3, 4)$, bude nabývat 
  hodnot menších než 6?    \
  4.  Předpokládejme, že $s_x^2=8,6$ při $15$ stupních volnosti. Doplňte následující
  výraz k určení 95\% intervalu spolehlivosti pro $\sigma_x^2$ (viz @eq-odhad-chisq-rozptyl).\        
  `c(15*` `r fitb(c("8.6", "8,6"))` `/qchisq(` `r fitb(c("0.975", "0,975"))`,  `r fitb(15)` `),` `15*` `r fitb(c("8.6", "8,6"))` `/qchisq(` `r fitb(c("0.025", "0,025"))`,  `15` `))`         
  5. Najděte následující hodnoty:
    - t_{0,95}(10), \chi_{0,1}^2(8), F_{0,01}(13;4)\         
:::