Cambiar a nuevo pipe

14-intro-bayesiana.Rmd

@@ -201,14 +201,14 @@ crear_verosim <- function(no_soles){
 # evaluar verosimilitud
 verosim <- crear_verosim(2)
 # ahora usamos regla de bayes para hacer tabla de probabilidades
-tabla_inferencia <- probs_inicial %>%
-  mutate(verosimilitud = map_dbl(theta, verosim)) %>%
-  mutate(inicial_x_verosim = prob_inicial * verosimilitud) %>%
+tabla_inferencia <- probs_inicial |>
+  mutate(verosimilitud = map_dbl(theta, verosim)) |>
+  mutate(inicial_x_verosim = prob_inicial * verosimilitud) |>
   # normalizar
   mutate(prob_posterior = inicial_x_verosim / sum(inicial_x_verosim))
 
-tabla_inferencia %>%
-  mutate(moneda_obs = moneda) %>%
+tabla_inferencia |>
+  mutate(moneda_obs = moneda) |>
   select(moneda_obs, theta, prob_inicial, verosimilitud, prob_posterior)
 ```
 
@@ -308,7 +308,7 @@ cualquier valor entre 0 y 1, pero es probable que tome un
 valor no tan lejano de 0.5. Por ejemplo, con probabilidad 0.95 creemos
 que $\theta$ está en el intervalo
 ```{r}
-quantile(sim_inicial$theta, c(0.025, 0.975)) %>% round(2)
+quantile(sim_inicial$theta, c(0.025, 0.975)) |> round(2)
 ```
 
 Es difícil justificar en abstracto por qué escogeriamos una inicial con esta
@@ -336,8 +336,8 @@ Concluimos entonces que la posterior tiene una distribución $\mathsf{Beta}(22,
 14)$. Podemos simular de la posterior usando código estándar para ver cómo luce:
 
 ```{r}
-sim_inicial <- sim_inicial %>% mutate(dist = "inicial")
-sim_posterior <- tibble(theta = rbeta(10000, 22, 14)) %>% mutate(dist = "posterior")
+sim_inicial <- sim_inicial |> mutate(dist = "inicial")
+sim_posterior <- tibble(theta = rbeta(10000, 22, 14)) |> mutate(dist = "posterior")
 sims <- bind_rows(sim_inicial, sim_posterior)
 ggplot(sims, aes(x = theta, fill = dist)) +
   geom_histogram(aes(x = theta), bins = 30, alpha = 0.5, position = "identity")
@@ -352,8 +352,8 @@ Podemos resumir de varias maneras. Por ejemplo, si queremos un estimador puntual
 usamos la media posterior:
 
 ```{r}
-sims %>% group_by(dist) %>%
-  summarise(theta_hat = mean(theta) %>% round(3))
+sims |> group_by(dist) |>
+  summarise(theta_hat = mean(theta) |> round(3))
 ```
 Nota que el estimador de máxima verosimilitud es $\hat{p} = 19/30 = 0.63$, que
 es ligeramente diferente de la media posterior. ¿Por qué?
@@ -363,8 +363,8 @@ suelen llamarse *intervalos de credibilidad*, por ejemplo:
 
 ```{r}
 f <- c(0.025, 0.975)
-sims %>% group_by(dist) %>%
-  summarise(cuantiles = quantile(theta, f) %>% round(2), f = f) %>%
+sims |> group_by(dist) |>
+  summarise(cuantiles = quantile(theta, f) |> round(2), f = f) |>
   pivot_wider(names_from = f, values_from = cuantiles)
 ```
 El segundo renglón nos da un intervalo posterior para $\theta$ de *credibilidad*
@@ -405,8 +405,8 @@ misma distribución inicial.
 puedes usar la aproximación de R, por ejemplo:
 
 ```{r}
-qbeta(0.025, shape1 = 22, shape2 = 14) %>% round(2)
-qbeta(0.975, shape1 = 22, shape2 = 14) %>% round(2)
+qbeta(0.025, shape1 = 22, shape2 = 14) |> round(2)
+qbeta(0.975, shape1 = 22, shape2 = 14) |> round(2)
 ```
 
 
@@ -435,14 +435,17 @@ $$P(\theta) = \frac{\alpha \theta_0^\alpha}{\theta^{\alpha + 1}}$$
 con soporte en $[\theta_0,\infty]$. Tenemos que escoger entonces el mínimo $\theta_0$ y
 el parámetro $\alpha$. En primer lugar, como sabemos que es una lotería nacional,
 creemos que no puede haber menos de unos 300 mil números, así que $\theta_0 = 300$.
-La función acumulada de la pareto es $1- (300/\theta)^\alpha$, así que el cuantil 99% es
+La función acumulada de la pareto es $1- (300/\theta)^\alpha$, así que si $\alpha = 1.1$
+el cuantil 99% es:
 
 ```{r}
 alpha <- 1.1
 (300/(0.01)^(1/alpha))
 ```
-es decir, alrededor de 20 millones de números.  Creemos que es un poco probable
+
+es decir, alrededor de 20 millones de números.  Creemos que es poco probable
 que el número de boletos sea mayor a esta cota.
+
 Nótese ahora que la posterior cumple (multiplicando verosimilitud por inicial):
 
 $$P(\theta|X_1,\ldots, X_n |\theta) \propto \theta^{-(n + 2.1)}$$
@@ -455,12 +458,12 @@ muestra de números:
 
 
 ```{r}
-loteria_tbl <- read_csv("data/nums_loteria_avion.csv", col_names = c("id", "numero")) %>%
+loteria_tbl <- read_csv("data/nums_loteria_avion.csv", col_names = c("id", "numero")) |>
   mutate(numero = as.integer(numero))
 set.seed(334)
-muestra_loteria <- sample_n(loteria_tbl, 25) %>%
+muestra_loteria <- sample_n(loteria_tbl, 25) |>
   mutate(numero = numero/1000)
-muestra_loteria %>% as.data.frame %>% head
+muestra_loteria |> as.data.frame() |> head()
 ```
 
 Podemos simular de una Pareto como sigue:
@@ -497,18 +500,19 @@ ggplot(sims_theta) +
   xlim(0, 15000) + scale_y_sqrt() +
   geom_rug(data = muestra_loteria, aes(x = numero))
 ```
+
 Nótese que cortamos algunos valores de la inicial en la cola derecha: un defecto
 de esta distribución inicial, con una cola tan larga a la derecha, es que
 pone cierto peso en valores que son poco creíbles y la vuelve poco apropiada para
-este problema. Regresamos más adelante a este problema.
+este problema. Regresaremos más adelante a este problema.
 
 Si obtenemos percentiles,
 obtenemos el intervalo
 
 ```{r}
 f <- c(0.025, 0.5, 0.975)
-sims_theta %>% group_by(dist) %>%
-  summarise(cuantiles = quantile(theta, f) %>% round(2), f = f) %>%
+sims_theta |> group_by(dist) |>
+  summarise(cuantiles = quantile(theta, f) |> round(2), f = f) |>
   pivot_wider(names_from = f, values_from = cuantiles)
 ```
 Estimamos entre 5.8 millones y 6.7 millones de boletos. El máximo en la muestra
@@ -643,6 +647,7 @@ no son tratables con análisis conjugado (veremos más adelante cómo tratarlos
 
 
 ### Ejemplo {-}
+
 Supongamos que queremos estimar la estatura de los cantantes de tesitura tenor con
 una muestra iid de tenores de Estados Unidos. Usaremos el modelo normal de forma que $X_i\sim \mathsf{N}(\mu, \sigma^2)$.
 
@@ -661,7 +666,7 @@ sigma_0 <- 7
 # disperisón
 a <- 3
 # ponemos 7 = sqrt(b/a) -> b = a * 64
-b <- a * sigma_0^2
+b <- a * sigma_0 ^ 2
 c(a = a, b = b)
 ```
 
@@ -707,7 +712,7 @@ la media poblacional de los cantantes.
 Podemos checar nuestros supuestos simulando posibles muestras usando
 sólo nuestra información previa:
 
-```{r, fig.width = 7, fig.height = 6}
+```{r, fig.width = 8, fig.height = 6}
 simular_normal_invgamma <- function(n, pars){
   mu_0 <- pars[1]
   n_0 <- pars[2]
@@ -721,8 +726,8 @@ simular_normal_invgamma <- function(n, pars){
   rnorm(n, mu, sigma)
 }
 set.seed(3461)
-sims_tbl <- tibble(rep = 1:20) %>%
-  mutate(estatura = map(rep, ~ simular_normal_invgamma(500, c(mu_0, n_0, a, b)))) %>%
+sims_tbl <- tibble(rep = 1:20) |>
+  mutate(estatura = map(rep, ~ simular_normal_invgamma(500, c(mu_0, n_0, a, b)))) |>
   unnest(cols = c(estatura))
 ggplot(sims_tbl, aes(x = estatura)) + geom_histogram() +
   facet_wrap(~ rep) +
@@ -741,9 +746,9 @@ Ahora podemos usar los datos para calcular nuestras posteriores.
 
 ```{r}
 set.seed(3413)
-cantantes <- lattice::singer %>%
-  mutate(estatura_cm = round(2.54 * height)) %>%
-  filter(str_detect(voice.part, "Tenor")) %>%
+cantantes <- lattice::singer |>
+  mutate(estatura_cm = round(2.54 * height)) |>
+  filter(str_detect(voice.part, "Tenor")) |>
   sample_n(20)
 cantantes
 ```
@@ -782,19 +787,19 @@ sim_params <- function(m, pars){
   a <- pars[3]
   b <- pars[4]
   # simular sigmas
-  sims <- tibble(tau = rgamma(m, a, b)) %>%
+  sims <- tibble(tau = rgamma(m, a, b)) |>
     mutate(sigma = 1 / sqrt(tau))
   # simular mu
-  sims <- sims %>% mutate(mu = rnorm(m, mu_0, sigma / sqrt(n_0)))
+  sims <- sims |> mutate(mu = rnorm(m, mu_0, sigma / sqrt(n_0)))
   sims
 }
-sims_inicial <- sim_params(5000, c(mu_0, n_0, a, b)) %>%
+sims_inicial <- sim_params(5000, c(mu_0, n_0, a, b)) |>
   mutate(dist = "inicial")
-sims_posterior <- sim_params(5000, pars_posterior) %>%
+sims_posterior <- sim_params(5000, pars_posterior) |>
   mutate(dist = "posterior")
 sims <- bind_rows(sims_inicial, sims_posterior)
 ggplot(sims, aes(x = mu, y = sigma, colour = dist)) +
-  geom_point()
+  geom_point(alpha = 0.4)
 ```
 
 Y vemos que nuestra posterior es consistente con la información inicial
@@ -803,32 +808,33 @@ ve como sigue. Hemos marcado también las medias posteriores de cada parámetro:
 media y desviación estándar.
 
 ```{r}
-medias_post <- sims %>% filter(dist == "posterior") %>%
-  select(-dist) %>%
+medias_post <- sims |> filter(dist == "posterior") |>
+  select(-dist) |>
   summarise(across(everything(), mean))
-ggplot(sims %>% filter(dist == "posterior"),
-    aes(x = mu, y = sigma)) +
+ggplot(sims |> filter(dist == "posterior"), aes(x = mu, y = sigma)) +
   geom_point(colour = "#00BFC4") +
   geom_point(data = medias_post, size = 5, colour = "black") +
   coord_equal()
 ```
+
 Podemos construir intervalos creíbles del 90% para estos dos parámetros, por ejemplo
 haciendo intervalos de percentiles:
 
 ```{r}
 f <- c(0.05, 0.5, 0.95)
-sims %>%
-  pivot_longer(cols = mu:sigma, names_to = "parametro") %>%
-  group_by(dist, parametro) %>%
-  summarise(cuantil = quantile(value, f) %>% round(1), f= f) %>%
+sims |>
+  pivot_longer(cols = mu:sigma, names_to = "parametro") |>
+  group_by(dist, parametro) |>
+  reframe(cuantil = quantile(value, f) |> 
+            round(1), f = f) |>
   pivot_wider(names_from = f, values_from = cuantil)
 ```
 
 Como comparación, los estimadores de máxima verosimlitud son
 
 ```{r}
 media_mv <- mean(cantantes$estatura_cm)
-sigma_mv <- mean((cantantes$estatura_cm - media_mv)^2) %>% sqrt
+sigma_mv <- mean((cantantes$estatura_cm - media_mv)^2) |> sqrt()
 c(media_mv, sigma_mv)
 ```
 
@@ -942,31 +948,33 @@ Y ahora simulamos otra muestra
 
 ```{r}
 muestra_sim <- simular_normal_invgamma(20, pars_posterior)
-muestra_sim %>% round(0)
+muestra_sim |> round(0)
 ```
 
 Podemos simular varias muestras y hacer una prueba de lineup:
 
 ```{r}
 library(nullabor)
 set.seed(9921)
-sims_obs <- tibble(.n = 1:19) %>%
-  mutate(estatura_cm = map(.n, ~ simular_normal_invgamma(20, pars_posterior))) %>%
+sims_obs <- tibble(.n = 1:19) |>
+  mutate(estatura_cm = map(.n, ~ simular_normal_invgamma(20, pars_posterior))) |>
   unnest(estatura_cm)
 pos <- sample(1:20, 1)
-lineup_tbl <- lineup(true = cantantes %>% select(estatura_cm),
+
+lineup_tbl <- lineup(true = cantantes |> select(estatura_cm),
                      samples = sims_obs, pos = pos)
 ggplot(lineup_tbl, aes(x = estatura_cm)) + geom_histogram(binwidth = 2.5) +
   facet_wrap(~.sample)
 ```
+
 Con este tipo de gráficas podemos checar desajustes potenciales de nuestro modelo.
 
 ```{block, type="ejercicio"}
 
 - ¿Puedes encontrar los datos verdaderos? ¿Cuántos seleccionaron los
 datos correctos?
-- Prueba hacer pruebas con una gráfica de cuantiles. ¿Qué problema
-ves y cómo lo resolverías?
+<!-- - Prueba hacer pruebas con una gráfica de cuantiles. ¿Qué problema -->
+<!-- ves y cómo lo resolverías? -->
 
 ```
 
@@ -997,8 +1005,8 @@ crear_sim_rep <- function(x){
   }
 }
 sim_rep <- crear_sim_rep(x)
-lineup_tbl <- map(1:5, ~ sim_rep(.x)) %>%
-  bind_rows() %>%
+lineup_tbl <- map(1:5, ~ sim_rep(.x)) |>
+  bind_rows() |>
   bind_rows(tibble(rep = 6, x_rep = x))
 ggplot(lineup_tbl, aes(x = x_rep)) +
   geom_histogram(bins = 15) +
@@ -1020,8 +1028,8 @@ crear_sim_rep <- function(x){
   }
 }
 sim_rep <- crear_sim_rep(x)
-lineup_tbl <- map(1:5, ~ sim_rep(.x)) %>%
-  bind_rows() %>%
+lineup_tbl <- map(1:5, ~ sim_rep(.x)) |>
+  bind_rows() |>
   bind_rows(tibble(rep = 6, x_rep = x))
 ggplot(lineup_tbl, aes(x = x_rep)) +
   geom_histogram(bins = 15) +
@@ -1051,13 +1059,13 @@ siempre, quisiéramos obtener un intervalo que exprese nuestra incertidumbre ace
 del valor que vamos a observar. Entonces haríamos:
 
 ```{r}
-sims_posterior <- sim_params(50000, pars_posterior) %>%
+sims_posterior <- sim_params(50000, pars_posterior) |>
   mutate(y_pred = rnorm(n(), mu, sigma))
-sims_posterior %>% head
+sims_posterior |> head()
 ```
 ```{r}
 f <- c(0.025, 0.5, 0.975)
-sims_posterior %>% summarise(f = f, y_pred = quantile(y_pred, f))
+sims_posterior |> summarise(f = f, y_pred = quantile(y_pred, f))
 ```
 
 Y con esto obtenemos el intervalo (163, 189), al 95%, para una nueva observación. Nótese
@@ -1105,7 +1113,7 @@ k_posterior <- nrow(muestra_loteria) + 1.1
 sims_pareto_posterior <- tibble(
   theta = rpareto(100000, lim_inf_post, k_posterior))
 # Simulamos una observación para cada una de las anteriores:
-sims_post_pred <- sims_pareto_posterior %>%
+sims_post_pred <- sims_pareto_posterior |>
   mutate(x_pred = map_dbl(theta, ~ runif(1, 0, .x)))
 # Graficamos
 ggplot(sims_post_pred, aes(x = x_pred)) +