---
title: "AMV1 Práctica 2.2 - Ej 8"
output:
  html_document:
    df_print: paged
---

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

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(magrittr)
library(glue)
library(stats)
library(tictoc)
library(gridExtra)
library(latex2exp)
library(tictoc)
library(gt)
library(ggforce)
```

# 8

Tenemos una muestra de t multivariadas:

$$
\mathbf{x}_i \overset{\text{iid}} \sim \mathcal{T}_{p,k} \\
\mathbf{x}_i = v^{1/2} \mathbf{z}_i \quad
\mathbf{z}_i \overset{\text{iid}} \sim \mathcal{N}_p(\mathbf{0}, \mathbf{I}_p) \quad
\frac{k}{v} \sim \chi^2_k
$$
$$
\mathbb{E}[\overline{\mathbf{x}}] = \mu \\
\mathbb{E}[\mathbf{S}] = \mathrm{Var}[\mathbf{x}] = \mathbb{E}v \cdot \Sigma =
\frac{k}{k - 2} \cdot \mathbf{I}_p \quad  k > 2
$$

Nos piden calcular, para cada muestra:

\begin{align*}

& D = || \overline{\mathbf{x}} − μ||^2
\overset{\mathbb{E}} \longrightarrow
\mathbf{0}
& k > 1
\\

& D_1 = \mathrm{Tr} ( \mathbf{S} \Sigma^{−1} )
\overset{\mathbb{E}} \longrightarrow
\frac{k}{k-2} \cdot p
& k > 2
\\

& D_2 = \frac{ \mathrm{Tr} ( \mathbf{S} \Sigma^{−1} ) }{ \mathbb{E}[v] }
\overset{\mathbb{E}} \longrightarrow
p
& k > 2

\end{align*}

### Mini Ejemplo

Para cada escenario $(k, n)$, tendremos una matriz $\mathbf{X} \in \mathbb{R}^{n \mathrm{x} 3}$, con los vectores $\mathbf{x}_i$ como filas.

Acá, ejemplo con $n = 10$ y $k = 2$:

```{r}
set.seed(12)

X <- mvtnorm::rmvt(10, df = 2, sigma = diag(3))
X
```

Calculamos el promedio de las $\mathbf{x}_i$ con `colMeans`:

```{r}
colMeans(X)
```

Calculamos $\mathbf{S}$ con `cov`:

```{r}
cov(X)
```

### Una muestra $\mathbf{x}_i \overset{\text{iid}} \sim \mathcal{T}_{3,k}$ ($i = 1, \dots, n$) para cada escenario $(n, k)$

```{r}
ESCENARIOS <-
  cross_df(list(
    n = c(10, 20, 50, seq(100, 1000, by = 50)),
    k = c(1, 2, 4, 10)
  ))

ESCENARIOS
```

```{r cache=T}
SIGMA <- diag(3)
MU <- rep(0, 3)

traza <- function(A) { sum(diag(A)) }

simular <- function() {
  ESCENARIOS %>%
    mutate(
      # X = map2(n, k, ~ mvtnorm::rmvt(n = ..1, df = ..2, sigma = SIGMA)),
      X = map2(n, k, mvtnorm::rmvt, sigma = SIGMA),
      x_barra = map(X, colMeans),
      S = map(X, cov),
      D = map_dbl(x_barra, ~ sum((. - MU)^2)),
      D_1 = map_dbl(S, ~ traza(. %*% solve(SIGMA))),
      # Esperanza_v = map_dbl(k, ~ ifelse(. > 2, ./(. - 2), NA)),
      Esperanza_v = ifelse(k > 2, k/(k - 2), NA),
      D_2 = D_1/Esperanza_v
    ) %>%
    select(-X, -x_barra, -S)
}

tic()
escenarios <- simular()
toc()

escenarios
```

```{r}
escenarios %<>% mutate_if(is.double, round, 3)

hacer_tabla_D <- function(.escenarios) {
  .escenarios %>%
    select(n, k, D) %>%
    spread(n, D) %>%
    gt() %>%
    tab_header(
      title = "D → 0; si k > 1"
    ) %>%
    tab_spanner(
      columns = matches("[0-9]"),
      label = "n"
    ) %>%
    fmt_number(
      columns = matches("[0-9]"),
      decimals = 3
    )
}

hacer_tabla_D1 <- function(.escenarios) {
  .escenarios %>%
    select(n, k, D_1) %>%
    spread(n, D_1) %>%
    gt() %>%
    tab_header(
      title = "D1 → (k / (k - 2)) p; si k > 2"
    ) %>%
    tab_spanner(
      columns = matches("[0-9]"),
      label = "n"
    ) %>%
    fmt_number(
      columns = matches("[0-9]"),
      decimals = 3
    )
}

hacer_tabla_D2 <- function(.escenarios) {
  .escenarios %>%
    select(n, k, D_2) %>%
    spread(n, D_2) %>%
    gt() %>%
    tab_header(
      title = "D2 → p; si k > 2"
    ) %>%
    tab_spanner(
      columns = matches("[0-9]"),
      label = "n"
    ) %>%
    fmt_number(
      columns = matches("[0-9]"),
      decimals = 3
    ) %>%
    fmt_missing(
      columns = everything(),
      missing_text = "-"
    )
}
```


```{r}
plot_facetado_k_n <- function(.escenarios) {
  .escenarios %>%
    mutate(k_str = fct_reorder(glue("k = {k}"), k)) %>%
    ggplot(aes(x = n)) +
    facet_grid(k_str ~ ., scales = "free_y")
}
```

```{r}
graficar_D <- function(.escenarios) {
  .escenarios %>%
    plot_facetado_k_n() +
    aes(y = D) +
    geom_hline(yintercept = 0, color = "Grey40", linetype = "dotted") +
    geom_line(color = "Grey20") +
    geom_point() +
    labs(
      title = TeX("Evolución de $D = || \\bar{x} - \\mu ||^2 \\rightarrow 0$"),
      subtitle = "Converge para k > 1"
    )
}
```

```{r}
graficar_D1 <- function(.escenarios) {
  plot_facetado_k_n(.escenarios) +
    aes(y = D_1) +
    geom_hline(
      mapping = aes(yintercept = Esperanza_v * 3),
      color = "ForestGreen", linetype = "dashed"
    ) +
    geom_point() +
    geom_line(color = "Grey50") +
    labs(
      title = TeX("Evolución de $D_1 = Tr ( S \\Sigma^{-1} ) \\rightarrow \\left( \\frac{k}{k - 2} \\right) \\cdot p$"),
      subtitle = "Converge si k > 2",
      y = TeX("$D_1$")
    )
}
```

```{r}
graficar_D2 <- function(.escenarios) {
  .escenarios %>%
    filter(!is.na(D_2)) %>%
    plot_facetado_k_n() +
    aes(y = D_2) +
    geom_hline(yintercept = 3, color = "ForestGreen", linetype = "dashed") +
    geom_point() +
    geom_line(color = "Grey50") +
    # geom_segment(aes(xend = n, yend = 3), color = "Red4") +
    labs(
      title = TeX("Evolución de $D_2 = Tr ( S \\Sigma^{-1} ) / E(v) \\rightarrow  p$"),
      subtitle = TeX("Sólo definido si k > 2"),
      y = TeX("$D_2$")
    )
}
```


```{r}
hacer_tabla_D(escenarios)
graficar_D(escenarios)
```

$$
\frac{k}{k - 2} \cdot p = \frac{4}{4 - 2} \cdot 3 = 6 \\
\frac{k}{k - 2} \cdot p = \frac{10}{10 - 2} \cdot 3 = 3.75
$$

```{r}
hacer_tabla_D1(escenarios)
graficar_D1(escenarios)
```


```{r}
hacer_tabla_D2(escenarios)
graficar_D2(escenarios)
```

## Consigna (b): $1.000$ iteraciones y promediar para cada $(n, k)$

```{r cache=T}
NITER <- 1000

tic()

simulaciones <-
  tibble(
    iter = 1:NITER,
    simulacion = map(iter, ~ simular())
  ) %>%
  unnest()

toc()
```

```{r message=FALSE, warning=FALSE}
promedios <-
  simulaciones %>%
  group_by(n, k, Esperanza_v) %>%
  summarise(
    count = n(),
    D = mean(D),
    D_1 = mean(D_1),
    D_2 = mean(D_2)
  ) %>%
  mutate_if(is.double, round, 3) %>%
  ungroup()

promedios %>% arrange(k, n)
```

```{r}
mensajito <- glue("Promedio de {NITER} iteraciones")

nota_al_pie <- function(data) {
  tab_source_note(
    data = data,
    source = mensajito
  )
}

caption <- labs(caption = mensajito)
```


```{r}
hacer_tabla_D(promedios) %>% nota_al_pie()
graficar_D(promedios) + caption
```


```{r fig.height=6, fig.width=6}
hacer_tabla_D1(promedios) %>% nota_al_pie()
graficar_D1(promedios) + caption
```


```{r}
hacer_tabla_D2(promedios) %>% nota_al_pie()
graficar_D2(promedios) + caption
```
