From c3ae06f2078f986b71afa9d47bea412c2dcaa821 Mon Sep 17 00:00:00 2001 From: tereom Date: Wed, 27 Nov 2024 15:20:28 -0600 Subject: [PATCH] Cambiar a nuevo pipe --- 14-intro-bayesiana.Rmd | 122 ++++++++++++++++++++++------------------- 1 file changed, 65 insertions(+), 57 deletions(-) diff --git a/14-intro-bayesiana.Rmd b/14-intro-bayesiana.Rmd index 1511e39..72160a5 100644 --- a/14-intro-bayesiana.Rmd +++ b/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,24 +808,25 @@ 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) ``` @@ -828,7 +834,7 @@ 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,7 +948,7 @@ 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: @@ -950,23 +956,25 @@ 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? + + ``` @@ -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)) +