06_odhad.qmd

# Bodový a intervalový odhad {#sec-odhad}

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

- Umět definovat vlastnosti odhadu.
- Naučit se aplikovat zákon velkých čísel.
- Provést bodový a intervalový odhad parametru.
- Užít $t$-test jako pomocný nástroj k sestavení odhadu.
:::

Obecně je odhadem v matematické statistice nazýváno určení parametru rozdělení hodnoty určitého znaku **základního souboru** s pomocí **výběrového souboru**.

Obecné charakteristiky základního souboru značíme písmeny řecké abecedy (např. $\mu$, $\sigma$, $\ldots$), pro výběrové charakteristiky volíme analogická písmena z latinky ($\bar{x}, s_x, \ldots$).

## Vlastnosti odhadu

V mnoha situacích potřebujeme odhadnout určitý parametr (střední hodnotu, $90\%$ kvantil atp.) neznámé náhodné veličiny. Tento parametr odhadujeme pomocí nějaké statistiky výběru z této veličiny. Například odhadujeme střední hodnotu veličiny pomocí průměru. Abychom mohli určit hodnotu parametru (např. střední hodnotu) přesně, musel by být výběr nekonečně velký. Jelikož toto v praxi nenastane, náš odhad je vždy více či méně odlišný od skutečné hodnoty parametru (např. střední hodnoty) neznámé veličiny. Důležitou vlastností odhadu je nestrannost - o nestranném odhadu mluvíme pokud střední hodnota odhadu je rovna neznámému parametru.

V přednáškách o charakteristikách náhodné veličiny a jejich odhadech je zmíněna rovnice

$$
\dfrac{\sum\limits_{i=1}^{n}(x_i - \bar{x})^2}{n-1}
$$

jako nestranný odhad veličiny $X$. Z definice rozptylu $\mathbb{E}[(X - \mathbb{E}(X^2))]$ však vyplývá odhad

$$
\dfrac{\sum\limits_{i=1}^{n}(x_i - \bar{x})^2}{n}.
$$

::: callout-tip
## Cvičení

  1. Doplňte následující kód.\
```{r, eval=FALSE}
X <- data.frame()
for (i in 1:nsim) {
  x <- rnorm(n)
  s_unb <- 1/(n-1)
  s_bia <- 1/n
  X <- rbind(X, c(s_unb, s_bia))
}
names(X) = c('UNB', 'BIA')
```
  a) Napište funkce pro odhad ropztylu.\
  b) Zamyslete se, které proměnné nejsou určeny a doplňte je.\
  c) Spočtěte pro každou metodu průměrný odhad a systematickou chybu tohoto odhadu.\
  d) Který odhad je nejméně vychýlený a v jaké situaci?\
:::


### Zákon velkých čísel

Zákon velkých čísel popisuje skutečnost, že s rostoucím počtem opakování nezávislých náhodných pokusů se empirické charakteristiky (realizované výběrové odhady), které popisují výsledky těchto pokusů blíží k teoretickým charakteristikám.

::: callout-tip
## Úloha

  1. Generujte s pomocí funkce `rnorm` postupně $10$, $10^2$, $10^3$ čísel se shodnou střední hodnotou a shodným rozptylem. Spočítejte $\bar{x}$ a ${s^2}$, použijte nápovědu.
Okomentujte výsledky. 
```{r, echo=FALSE, fig.align='center', fig.width=10}
par(mfrow = c(1, 3))
set.seed(13)
x10 <- rnorm(10, 0, 1)
set.seed(13)
x100 <- rnorm(100, 0, 1)
set.seed(13)
x1000 <- rnorm(1000, 0, 1)
curve(dnorm(x, mean = 0 , sd = 1), from = -4, to = 4, xlab = "", ylab = "", main = "A)", cex.main = 2)
rug(x10, lwd = 2, col = "#00000025")
rug(mean(x10), lwd = 2, col = "orangered")
segments(x0 = mean(x10), y0 = 0, x1 = mean(x10), y1 = dnorm(mean(x10), 0, 1), col = "orangered")
segments(x0 = 0, y0 = 0, x1 = 0, y1 = dnorm(0, 0, 1), lty = "dashed", col = "gray25")
curve(dnorm(x, mean = 0 , sd = 1), from = -4, to = 4, xlab = "", ylab = "", main = "B)", cex.main = 2)
rug(x100, lwd = 2, col = "#00000025")
rug(mean(x100), lwd = 2, col = "orangered")
segments(x0 = mean(x100), y0 = 0, x1 = mean(x100), y1 = dnorm(mean(x100), 0, 1), col = "orangered")
segments(x0 = 0, y0 = 0, x1 = 0, y1 = dnorm(0, 0, 1), lty = "dashed", col = "gray25")
curve(dnorm(x, mean = 0 , sd = 1), from = -4, to = 4, xlab = "", ylab = "", main = "C)", cex.main = 2)
rug(x1000, lwd = 2, col = "#00000025")
rug(mean(x1000), lwd = 2, col = "orangered")
segments(x0 = mean(x1000), y0 = 0, x1 = mean(x1000), y1 = dnorm(mean(x1000), 0, 1), col = "orangered")
segments(x0 = 0, y0 = 0, x1 = 0, y1 = dnorm(0, 0, 1), lty = "dashed", col = "gray25")
```
:::

## Bodový odhad

Bodovým odhadem se se rozumí jednočíselná hodnota, která reprezentuje vybraný moment statistického souboru jako celek. Bodovým odhadem je například výběrový průměr nebo výběrový rozptyl.

```{r, collapse=TRUE}
set.seed(100)
x <- rnorm(10, 10, 10)
mean(x)
var(x)
```


### Odhad parametru $\mu$, neboli střední hodnoty normálního rozdělení s neznámým rozptylem

Takto formulovaný bodový a intervalový je jednou z nejčastěji prováděných úloh. Nejprve bodovým odhadem zjistíme **výběrový průměr** souboru.

```{r}
mean(x)
```

Vidíme, že $x =$ `r mean(x)`. Pro tento průmer následně spočítáme interval spolehlivosti.

## Intervalový odhad

V případech, kdy chceme znát polohu bodového odhadu s nějakou danou pravděpodobností, můžeme se pokusit zkonstruovat tzv. **intervalový odhad**. @eq-odhad-intervalovyodhad je rovnice pro intervalový odhad střední hodnoty při neznámém rozptylu.


$$ 
\bar{X} \pm \dfrac{s}{\sqrt{n}}t_{1-\alpha/2}(n-1)
$$ {#eq-odhad-intervalovyodhad}

$100(1-\alpha)\%$ interval spolehlivosti je rozmezí, ve kterém se usuzovaná hodnota základního souboru bude nacházet s určitou pravěděpodobností. Je tomu tak současně za předpokladu, že samotná poloha
hledaného parametru je konstantní.

```{r}
alpha <- 0.05
cbind(
  mean(x) - sd(x)/length(x)*qt(p = 1 - alpha/2, df = length(x) - 1),
  mean(x) + sd(x)/length(x)*qt(p = 1 - alpha/2, df = length(x) - 1)
)
```

## Maximálně věrohodný odhad

Maximálně věrohodný odhad je takový, který maximalizuje věrohodnostní funkci $L$ pro výběr $(x_1, x_2, \ldots, x_n)$

$$
L(\theta, x_1, x_2, \ldots, x_n) = p(x_1, \theta) \cdot p(x_2, \theta) \cdot \ldots p(x_n, \theta) = \prod\limits_{i=1}^{n}p(x_i, \theta)
$$ Konstrukce věrohodnostní funkce

### Maximálně věrohodný odhad pro normální rozdělení

Vjděme z pravděpodobnostní funkce normálního rozdělení s předpisem @eq-maxlike-normalni-rozdeleni, který 

$$
f(y|\mu, \sigma) = \dfrac{1}{\sigma\sqrt{2\pi}}\exp\left[-\dfrac{(y-\mu)^2}{2\sigma^2}\right]
$$ {#eq-maxlike-normalni-rozdeleni}

Věrohodnostní funkce je vyjádřena:

$$
L\left(\mu, \sigma^2; x_1, \ldots, x_n\right) = (2\pi\sigma^2)^{-n/2}\exp\left(-\dfrac{1}{2\sigma^2}\sum_{j=1}^{n}(x_j-\mu)^2\right)
$$ {#eq-05odhad-likelihood}

Dva parametry $\mu$ a $\sigma$ jsou neznámé

```{r}
mu <- 5 # <1>
sigma <- 1 # <1>

n <- 1000 # <2>
x <- rnorm(n, mu, sigma^2) # <2>

norm.lik <- function(pars, x) { # <3>
  n <- length(x) # <3>
  mu <- pars[1]
  sigma2 <- pars[2]
  logl <- -0.5 * n * log(2*pi) - 0.5 * n * log(sigma2) - 
    (1/(2 * sigma2))*sum((x - mu)^2) # <3>
  return(-logl)# <4>
}


result <- optim(par = c(0.1, 2), #<5>
                fn = norm.lik, #<5>
                x = x, #<5>
                method = "BFGS", #<5>
                hessian = TRUE) #<5>

mle_mean <- result$par[1] #<6>
mle_sigma2 <- result$par[2] #<6>
```
1. Skutečné parametry $\mu$ a $\sigma$, které použijeme ke generování náhodných čísel.
2. Získáme pseudonáhodná čísla z Normálního rozdělení.
3. Definice věrohodnostní funkce @eq-05odhad-likelihood,
4. je nutné pro optimalizační algoritmus vrátit v záporné hodnotě.
5. S použitím funkce `optim()` hledáme maximálně věrohodný odhad. Je potřeba zadat
počáteční parametry, optimalizovanou funkci, vektor měření a methodu řešení. Pokud 
bychom chtěli i intervalový odhad, je nutné nechat spočítat hessián, neboli determinant
Hessovy matice.
6. Uložíme získané hodnoty odhadů.

MLE pro střední hodnotu: `r mle_mean`\     
MLE pro rozptyl: `r mle_sigma2`.

## $t$-test {#ttest}
Testování hypotéz podrobněji probereme v @sec-testovani. Nyní si nicméně ukážeme, 
že výstupem funkce `t.test()` je rovněž intervalový odhad pro střední hodnotu Normálního
rozdělení při neznámém rozptylu.
Testovací statistika pro oboustrannou alternativu má hodnotu

$$
\dfrac{|\bar{x} - \mu_0|}{s}\sqrt{n} > t_{\alpha/2}(n-1)
$$ a pro jednostrannou alternativu $\mu > \mu_0$

$$
\dfrac{\bar{x} - \mu_0}{s}\sqrt{n} > t_{\alpha}(n-1)
$$ respektive

$$
\dfrac{\bar{x} - \mu_0}{s}\sqrt{n} < t_{\alpha}(n-1)
$$ pro $\mu < \mu_0$. $n-1$ je počet stupňů volnosti.

::: callout-tip
## Úloha

  1. Spočítejte pomocí funkce `t.test` intervalový odhad pro `x = rnorm(100)` a `set.seed(100)`.\
:::