From fafe41a8a96de37d9d05205bb0b642525b1199ec Mon Sep 17 00:00:00 2001 From: tereom Date: Thu, 12 Dec 2024 11:09:55 -0600 Subject: [PATCH] Agregar notas stan --- 16-bayes-mcmc.Rmd | 1 + 17-inferencia-stan-01.Rmd | 422 +++++++++++++++++++++++++++++++ stan/modelo-1.stan | 20 ++ stan/modelo-2.stan | 20 ++ stan/modelo-3.stan | 20 ++ stan/modelo-cantantes.stan | 25 ++ stan/modelo-examenes-no-inf.stan | 24 ++ stan/modelo-examenes.stan | 23 ++ 8 files changed, 555 insertions(+) create mode 100644 17-inferencia-stan-01.Rmd create mode 100644 stan/modelo-1.stan create mode 100644 stan/modelo-2.stan create mode 100644 stan/modelo-3.stan create mode 100644 stan/modelo-cantantes.stan create mode 100644 stan/modelo-examenes-no-inf.stan create mode 100644 stan/modelo-examenes.stan diff --git a/16-bayes-mcmc.Rmd b/16-bayes-mcmc.Rmd index 9ade9b6..f66c51e 100644 --- a/16-bayes-mcmc.Rmd +++ b/16-bayes-mcmc.Rmd @@ -746,6 +746,7 @@ g_sim <- ggplot(sims_tbl %>% filter(iter_num < 3000), aes(x = iter_num, y = theta)) + geom_line() + ylim(c(0, 0.5)) g_sim + dist_bplot + plot_layout(widths = c(5, 1)) ``` + Donde vemos que esta cadena parece mezclar bien (está explorando la totalidad de la distribución objetivo), y también parece estar en un estado estable. diff --git a/17-inferencia-stan-01.Rmd b/17-inferencia-stan-01.Rmd new file mode 100644 index 0000000..e380b0a --- /dev/null +++ b/17-inferencia-stan-01.Rmd @@ -0,0 +1,422 @@ +# Ejemplos de inferencia bayesiana en Stan I + +```{r, echo=FALSE} +knitr::opts_chunk$set(error = TRUE) +``` + + +En esta parte veremos cómo correr y diagnosticar en Stan varios ejemplos vistos en clase. Para instalar *cmdstanr* y Stan, puedes ver [aquí](https://mc-stan.org/cmdstanr/). En python puedes usar *pystan*, por ejemplo. + +```{r, message=FALSE} +# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) +library(cmdstanr) +library(posterior) +library(tidyverse) +``` + + +## Estimación de una proporción {-} + +Escribimos el código para el modelo en un archivo *modelo-1.stan*, +y compilamos: + +```{r, message=FALSE} +archivo_stan <- file.path("stan/modelo-1.stan") +# compilar +mod <- cmdstan_model(archivo_stan) +``` + + +```{r, message=FALSE} +mod +``` + +Pasamos datos y muestreamos + +```{r} +datos_lista <- list(n = 30, y = 19) +ajuste <- mod$sample( + data = datos_lista, + seed = 1234, + chains = 4, + parallel_chains = 4, + refresh = 500) +``` + +Checamos diagnósticos: + +```{r} +ajuste$cmdstan_diagnose() +``` + +Si no hay problemas, podemos ver el resumen: + +```{r} +ajuste$summary() +``` + +Donde verificamos que el tamaño de muestra efectivo (ess) y el diagnóstico de +$\hat{R}$ son apropiados. + +Podemos ver las cadenas de la siguiente forma: + +```{r} +theta_tbl <- ajuste$draws(c("theta", "theta_inicial")) %>% as_draws_df() +ggplot(theta_tbl, aes(x = .iteration, y = theta)) + + geom_line() + + facet_wrap(~.chain, ncol = 1) +``` + +Y replicamos la gráfica de las notas haciendo: + +```{r, fig.width=6, fig.height = 4} +sims_tbl <- theta_tbl %>% pivot_longer(theta:theta_inicial, names_to = "dist", values_to = "theta") +ggplot(sims_tbl, aes(x = theta, fill = dist)) + + geom_histogram(aes(x = theta), bins = 30, alpha = 0.5, position = "identity") +``` + +## Estimación del máximo de una uniforme {-} + +Tomamos el ejemplo de los boletos de lotería, + +```{r} +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) %>% + mutate(numero = numero/1000) +``` + +```{r, message=FALSE, include=FALSE} +archivo_stan <- file.path("stan/modelo-2.stan") +# compilar +mod <- cmdstan_model(archivo_stan) +``` + + +```{r, message=FALSE} +mod +``` + +Pasamos datos y muestreamos: + +```{r} +datos_lista <- list(n = nrow(muestra_loteria), y = muestra_loteria$numero) +ajuste <- mod$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_warmup = 5000, + iter_sampling = 20000, + parallel_chains = 4, + refresh = 5000) +``` + +Checamos diagnósticos: + +```{r} +ajuste$cmdstan_diagnose() +``` + +Si no hay problemas, podemos ver el resumen: + +```{r} +resumen <- ajuste$summary() +resumen +``` + +El intervalo 95% que obtenemos es: + +```{r} +ajuste$draws("theta") %>% as_draws_df() %>% + summarise(theta_inf = quantile(theta, 0.025), + theta_mediana = quantile(theta, 0.5), + theta_sup = quantile(theta, 0.975)) +``` + +Podemos ahora intentar con la inicial gamma que nos pareció +más intuitiva, aún cuando el modelo no es conjugado: + +```{r, message=FALSE} +archivo_stan <- file.path("stan/modelo-3.stan") +# compilar +mod <- cmdstan_model(archivo_stan) +``` + + +```{r, message=FALSE} +mod +``` + +Pasamos datos y muestreamos + +```{r} +datos_lista <- list(n = nrow(muestra_loteria), y = muestra_loteria$numero) +ajuste <- mod$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_sampling = 10000, + parallel_chains = 4, + refresh = 2000) +``` + +Checamos diagnósticos: + +```{r} +ajuste$cmdstan_diagnose() +``` + +Si no hay problemas, podemos ver el resumen: + +```{r} +resumen <- ajuste$summary() +resumen +``` + +El intervalo 95% que obtenemos es: + +```{r} +ajuste$draws("theta") %>% as_draws_df() %>% + summarise(theta_inf = quantile(theta, 0.025), + theta_mediana = quantile(theta, 0.5), + theta_sup = quantile(theta, 0.975)) +``` + +Y la posterior se ve como sigue: + +```{r} +theta_post_sim <- ajuste$draws("theta") %>% as.numeric +qplot(theta_post_sim) +``` + +## Ejemplo de cantantes {-} + +Haremos el ejemplo no conjugado de estaturas de cantantes, con +verificación predictiva posterior. + +```{r, message=FALSE, include=FALSE} +archivo_stan <- file.path("stan/modelo-cantantes.stan") +# compilar +mod <- cmdstan_model(archivo_stan) +``` + + +```{r, message=FALSE} +mod +``` + +Pasamos datos y muestreamos + +```{r} +set.seed(3413) +cantantes <- lattice::singer %>% + mutate(estatura_cm = (2.54 * height)) %>% + filter(str_detect(voice.part, "Tenor")) %>% + sample_n(20) +datos_lista <- list(n = nrow(cantantes), y = cantantes$estatura_cm) +ajuste <- mod$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_warmup = 4000, + iter_sampling = 4000, + refresh = 2000) +``` + +Checamos diagnósticos: + +```{r} +ajuste$cmdstan_diagnose() +``` + +Si no hay problemas, podemos ver el resumen: + +```{r} +resumen <- ajuste$summary() +resumen +``` + +El intervalo 95% que obtenemos es: + +```{r} +ajuste$draws(c("mu", "sigma")) %>% as_draws_df() %>% +ggplot(aes(x = mu, y = sigma)) + geom_point(alpha = 0.1) + + coord_equal() +``` + + +Y ahora extraemos algunas replicaciones de la posterior predictiva: + +```{r} +y_sim_tbl <- ajuste$draws("y_sim") %>% as_draws_df() %>% + pivot_longer(cols = starts_with("y_sim"), "nombre") %>% + separate(nombre, c("nombre", "n_obs", "vacio"), "[\\[\\]]") %>% + select(-nombre, -vacio) %>% + filter(.chain == 1, .iteration < 12) %>% + select(.iteration, value) %>% + bind_rows(tibble(.iteration = 12, value = round(cantantes$estatura_cm, 0))) +ggplot(y_sim_tbl, aes(sample = value)) + + geom_qq() + + facet_wrap(~ .iteration) +``` + + +### Ejemplo: exámenes + +```{r} +sim_formas <- function(p_azar, p_corr){ + tipo <- rbinom(1, 1, 1 - p_azar) + if(tipo==0){ + # al azar + x <- rbinom(1, 10, 1/5) + } else { + # no al azar + x <- rbinom(1, 10, p_corr) + } + x +} +set.seed(12) +muestra_1 <- map_dbl(1:200, ~ sim_formas(0.35, 0.5)) +``` + + + +```{r, message = FALSE, include=FALSE} +archivo_stan <- file.path("stan/modelo-examenes.stan") +# compilar +mod_informado <- cmdstan_model(archivo_stan) +``` + + +```{r, message = FALSE} +mod_informado +``` + +Pasamos datos y muestreamos + +```{r} +set.seed(3413) +datos_lista <- list(n = length(muestra_1), y = muestra_1) +ajuste <- mod_informado$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_warmup = 4000, + iter_sampling = 4000, + refresh = 2000) +``` +Checamos diagnósticos: + +```{r} +ajuste$cmdstan_diagnose() +``` + +Si no hay problemas, podemos ver el resumen: + +```{r} +resumen <- ajuste$summary() +resumen +``` + +```{r} +sims_theta_tbl <- + ajuste$draws(c("theta_azar", "theta_corr")) %>% + as_draws_df() +ggplot(sims_theta_tbl, aes(x = theta_azar, y = theta_corr)) + + geom_point(alpha = 0.1) +``` + + + +## Ejemplo: exámenes, mal identificado + +En este ejemplo cambiamos los parámetros de la simulación +y ponemos iniciales poco informativas. + +```{r} +set.seed(12) +muestra_2 <- map_dbl(1:200, ~ sim_formas(0.35, 0.21)) +``` + + + +```{r, message = FALSE, include=FALSE} +archivo_stan <- file.path("stan/modelo-examenes-no-inf.stan") +# compilar +mod_no_inf <- cmdstan_model(archivo_stan) +``` + + +```{r, message = FALSE} +mod_no_inf +``` + +Pasamos datos y muestreamos + +```{r} +set.seed(3413) +datos_lista <- list(n = length(muestra_2), y = muestra_2) +ajuste <- mod_no_inf$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_warmup = 4000, + iter_sampling = 4000, + refresh = 2000) +``` +Y obtenemos mensajes de advertencia (divergencias), lo que implica que +los diagnósticos indican que es posible que las cadenas no hayan explorado +suficientemente la posterior + +```{r} +ajuste$cmdstan_diagnose() +``` + + +```{r} +sims_theta_tbl <- + ajuste$draws(c("theta_azar", "theta_corr")) %>% + as_draws_df() +ggplot(sims_theta_tbl, aes(x = theta_azar, y = theta_corr)) + + geom_point(alpha = 0.1) +``` + +Donde vemos que el problema es serio: cuando $\theta_{azar}$ es chico, los +datos son consistentes con valores de $\theta_{corr}$ cercanos a 0.2. Pero +también es posible que $\theta_{azar}$ sea grande, y en ese caso tenemos +poca información acerca de $\theta_corr$. + +- La geometría de esta posterior +hace difícil para el algoritmo establecer el tamaño de paso correcto para explorar +adecuadamente esta posterior (con un tamaño de paso chico, explora lentamente, pero +si se hace más grande, entonces aparecen problemas numéricos y rechazos en los +"embudos" de esta posterior). +- Sin embargo, este resultado generalmente lo rechazaríamos, pues sabemos que +una proporción considerable de estudiantes **no** contesta al azar. + +Si corremos el modelo informado con la muestra de esta población, obtenemos: + + + +```{r} +set.seed(3413) +datos_lista <- list(n = length(muestra_2), y = muestra_2) +ajuste <- mod_informado$sample( + data = datos_lista, + seed = 1234, + chains = 4, + iter_warmup = 4000, + iter_sampling = 4000, + refresh = 2000) +``` + +No tenemos problemas numéricos, y la posterior se ve como sigue: + +```{r} +sims_theta_tbl <- + ajuste$draws(c("theta_azar", "theta_corr")) %>% + as_draws_df() +ggplot(sims_theta_tbl, aes(x = theta_azar, y = theta_corr)) + + geom_point(alpha = 0.1) +``` diff --git a/stan/modelo-1.stan b/stan/modelo-1.stan new file mode 100644 index 0000000..30c84b0 --- /dev/null +++ b/stan/modelo-1.stan @@ -0,0 +1,20 @@ +// Ejemplo de estimación de una proporcion +data { + int n; // número de pruebas + int y; //numero de éxitos y fracasos +} + +parameters { + real theta; +} + +model { + // inicial + theta ~ beta(3, 3); + y ~ binomial(n, theta); +} + +generated quantities { + real theta_inicial; + theta_inicial = beta_rng(3, 3); +} diff --git a/stan/modelo-2.stan b/stan/modelo-2.stan new file mode 100644 index 0000000..a5a86a4 --- /dev/null +++ b/stan/modelo-2.stan @@ -0,0 +1,20 @@ +// Ejemplo de estimación del máximo de uniforme +data { + int n; // número de observaciones + array[n] real y; //datos observados +} + +transformed data{ + real y_max; + y_max = max(y); +} +parameters { + real theta; +} + +model { + // inicial + theta ~ pareto(300, 1.1); + y ~ uniform(0, theta); +} + diff --git a/stan/modelo-3.stan b/stan/modelo-3.stan new file mode 100644 index 0000000..8a00a01 --- /dev/null +++ b/stan/modelo-3.stan @@ -0,0 +1,20 @@ +// Ejemplo de estimación del máximo de uniforme +data { + int n; // número de observaciones + array[n] real y; //datos observados +} + +transformed data{ + real y_max; + y_max = max(y); +} +parameters { + real theta; +} + +model { + // inicial + theta ~ gamma(5, 0.0001); + y ~ uniform(0, theta); +} + diff --git a/stan/modelo-cantantes.stan b/stan/modelo-cantantes.stan new file mode 100644 index 0000000..6e6e986 --- /dev/null +++ b/stan/modelo-cantantes.stan @@ -0,0 +1,25 @@ +// Ejemplo de modelo normal para estaturas de cantantes +data { + int n; // número de observaciones + array[n] real y; //datos observados +} + +parameters { + real mu; + real sigma; +} + +model { + // inicial + mu ~ normal(175, 3); + sigma ~ uniform(2, 20); + y ~ normal(mu, sigma); +} + +generated quantities { + array[n] real y_sim; + for(i in 1:n){ + y_sim[i] = normal_rng(mu, sigma); + } + +} diff --git a/stan/modelo-examenes-no-inf.stan b/stan/modelo-examenes-no-inf.stan new file mode 100644 index 0000000..cb1d79d --- /dev/null +++ b/stan/modelo-examenes-no-inf.stan @@ -0,0 +1,24 @@ +// Ejemplo de estimación del máximo de uniforme +data { + int n; // número de observaciones + array[n] int y; //datos observados +} + + +parameters { + real theta_azar; + real theta_corr; +} + +model { + // inicial + theta_azar ~ beta(1, 1); + theta_corr ~ beta(1, 1); + // en este caso, agregamos términos directamente a la log posterior + for(i in 1:n){ + target+= log_sum_exp( + log(theta_azar) + binomial_lpmf(y[i] | 10, 0.20), + log(1 - theta_azar) + binomial_lpmf(y[i] | 10, theta_corr)); + } +} + diff --git a/stan/modelo-examenes.stan b/stan/modelo-examenes.stan new file mode 100644 index 0000000..9332998 --- /dev/null +++ b/stan/modelo-examenes.stan @@ -0,0 +1,23 @@ +// Ejemplo de estimación del máximo de uniforme +data { + int n; // número de observaciones + array[n] int y; //datos observados +} + + +parameters { + real theta_azar; + real theta_corr; +} + +model { + // inicial + theta_azar ~ beta(1, 5); + theta_corr ~ beta(7, 3); + // en este caso, agregamos términos directamente a la log posterior + for(i in 1:n){ + target+= log_sum_exp( + log(theta_azar) + binomial_lpmf(y[i] | 10, 0.20), + log(1 - theta_azar) + binomial_lpmf(y[i] | 10, theta_corr)); + } +}