From 5460c2ae73e36165ed0aec531a34c9b1cb4cb4ae Mon Sep 17 00:00:00 2001 From: tereom Date: Wed, 27 Nov 2024 15:53:51 -0600 Subject: [PATCH] Agregar notas calibracion --- 15-bayesiana-calibracion.Rmd | 870 +++++++++++++++++ 16-bayes-mcmc.Rmd | 1771 ++++++++++++++++++++++++++++++++++ renv.lock | 11 + 3 files changed, 2652 insertions(+) create mode 100644 15-bayesiana-calibracion.Rmd create mode 100644 16-bayes-mcmc.Rmd diff --git a/15-bayesiana-calibracion.Rmd b/15-bayesiana-calibracion.Rmd new file mode 100644 index 0000000..4e94138 --- /dev/null +++ b/15-bayesiana-calibracion.Rmd @@ -0,0 +1,870 @@ +# Calibración bayesiana y Regularización + +```{r setup, include=FALSE, message=FALSE} +library(tidyverse) +library(patchwork) +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center', fig.width = 6, fig.height=4) +comma <- function(x) format(x, digits = 2, big.mark = ",") +theme_set(theme_minimal()) +``` + +El enfoque bayesiano se puede formalizar coherentemente en términos +de probabilidades subjetivas, y como vimos, esta es una fortaleza del enfoque bayesiano. + +En la práctica, sin embargo, muchas veces puede ser difícil argumentar en términos +exclusivos de probabilidad subjetiva, aunque hagamos los esfuerzos apropiados para +incorporar la totalidad de información que distintos actores involucrados pueden tener. + +Consideremos, por ejemplo, que INEGI produjera un intervalo creíble del 95% +para el ingreso mediano de los hogares de México. Aún cuando nuestra metodología sea +transparente y correctamente informada, algunos investigadores interesados +puede ser que tengan recelo en usar esta información, y quizá preferirían +hacer estimaciones propias. Esto restaría valor al trabajo cuidadoso que pusimos +en nuestras estimaciones oficiales. + +Por otra parte, el enfoque frecuentista provee de ciertas garantías mínimas +para la utilización de las estimaciones, que no dependen de la interpretación +subjetiva de la probabilidad, sino de las propiedades del muestreo. +Consideremos la cobertura de los intervalos de confianza: + +- Bajo ciertos supuestos de nuestros modelos, la probabilidad de que +un intervalo de confianza del 95% cubra al verdadero valor poblacional es del 95%. Esta +probabilidad es sobre las distintas muestras que se pueden obtener según el +diseño del muestreo. + +Los intervalos creíbles en principio no tienen por qué cumplir esta propiedad, pero +consideramos que en la práctica es una garantía mínima que deberían cumplir. + +El enfoque resultante se llama **bayesiano calibrado**, @little2011 . La idea +es seguir el enfoque bayesiano usual para construir nuestras estimaciones, +pero verificar hasta donde sea posible que los intervalos resultantes satisfacen +alguna garantía frecuentista básica. + +**Observación**. checar que la cobertura real es similar a la nominal es importante +en los dos enfoques: frecuentista y bayesiano. Los intervalos frecuentistas, como +hemos visto, generalmente son aproximados, y por lo tanto no cumplen automáticamente +esta propiedad de calibración. + +## Enfoque bayesiano y frecuentista {-} + +Los métodos estadísticos clásicos toman el punto de vista **frecuentista** y se +basa en los siguientes puntos (@Wasserman): + +1. La probabilidad se interpreta como un límite de frecuencias relativas, donde +las probabilidades son propiedades objetivas en el mundo real. + +2. En un modelo, los parámetros son constantes fijas (desconocidas). Como +consecuencia, no se pueden realizar afirmaciones probabilísticas útiles en +relación a éstos. + +3. Los procedimientos estadísticos deben diseñarse con el objetivo de tener +propiedades frecuentistas bien definidas. Por ejemplo, un intervalo de confianza +del $95\%$ debe contener el verdadero valor del parámetro con frecuencia límite +de al menos el $95\%$. + +En contraste, el acercamiento **Bayesiano** muchas veces se describe por los +siguientes postulados: + +1. La probabilidad describe grados de creencia, no frecuencias limite. Como +tal uno puede hacer afirmaciones probabilísticas acerca de muchas cosas y no +solo datos sujetos a variabilidad aleatoria. Por ejemplo, puedo decir: "La +probabilidad de que Einstein tomara una taza de té el primero de agosto de $1948$" +es $0.35$, esto no hace referencia a ninguna frecuencia relativa sino que refleja +la certeza que yo tengo de que la proposición sea verdadera. + +2. Podemos hacer afirmaciones probabilísticas de parámetros. + +3. Podemos hacer inferencia de un parámetro $\theta$ por medio de +distribuciones de probabilidad. Las inferencias como estimaciones puntuales y +estimaciones de intervalos se pueden extraer de dicha distribución. + +Finalmente, en el enfoque **bayesiano calibrado** (@little2011): + +1. Usamos el enfoque bayesiano para modelar y hacer afirmaciones probabilísticas +de los parámetros. +2. Buscamos cumplir las garantías frecuentistas del inciso 3). + +## Ejemplo: estimación de una proporción {-} + +Recordamos nuestro problema de estimación de una proporcion $\theta$. Usando +la distribución inicial $p(\theta)\sim \mathsf{Beta}(2,2)$, y la verosimilitud estándar +binomial, vimos que la posterior cuando observamos $k$ +éxitos es $$p(\theta|k) \sim \mathsf{Beta}(k + 2, n - k + 2)$$. + +La media posterior es +$$\frac{k + 2}{n + 4} $$ +que podemos interpretar como: agrega 2 éxitos y 2 fracasos a los datos observados +y calcula la proporción de éxitos. Un intervalo posterior de credibilidad del 95% se calcula +encontrando los cuantiles 0.025 y 0.975 de una $\mathsf{Beta}(k + 2, n - k + 2)$ + +$$I_a = \left [q_{0.025}(k+2, n+4), q_{0.975}(k+2, n+4)\right ]$$ +Que compararemos con el intervalo usual de Wald: si $\hat{\theta} = \frac{k}{n}$, entonces + +$$I_w = \left [\hat{\theta} - 2 \sqrt{\frac{\hat{\theta}(1-\hat{\theta})}{n}}, \hat{\theta} + 2 \sqrt{\frac{\hat{\theta}(1-\hat{\theta})}{n}}\right]$$ +¿Cómo podemos comparar la calibración de estos dos intervalos? Nominalmente, deben +tener cobertura de 95%. Hagamos un ejercicio de simulación para distintos +tamaños de muestra $n$ y posibles valores $\theta\in (0,1)$: + +```{r} +set.seed(332) +simular_muestras <- function(M, n, p){ + k = rbinom(M, n, p) + tibble(rep = 1:M, n = n, p = p, k = k) +} +intervalo_wald <- function(n, k){ + p_hat <- k / n + ee_hat <- sqrt(p_hat * (1 - p_hat) / n) + tibble(inf = p_hat - 2 * ee_hat, sup = p_hat + 2 * ee_hat) +} +intervalo_bayes <- function(n, k, a = 2, b = 2){ + a <- k + a + b <- n - k + b + tibble(inf = qbeta(0.025, a, b), sup = qbeta(0.975, a, b)) +} +set.seed(812) +ejemplo <- simular_muestras(5, 20, 0.4) +``` + +```{r} +ejemplo |> mutate(intervalo = intervalo_wald(n, k)) |> pull(intervalo) |> + bind_cols(ejemplo) |> select(-rep) +``` + +```{r} +ejemplo |> mutate(intervalo = intervalo_bayes(n, k)) |> pull(intervalo) |> + bind_cols(ejemplo) |> select(-rep) +``` + +¿Cuáles de estos intervalos cubren al verdadero valor? Nótese que **no podemos +descalificar a ningún método por no cubrir una vez**. Es fácil producir un intervalo +con 100% de cobertura: (0,1). Pero no nos informa dónde es probable que esté +el parámetro. + +Sin embargo, podemos checar la cobertura frecuentista haciendo una cantidad grande de simulaciones: + +```{r} +parametros <- crossing(n = c(5, 10, 30, 60, 100, 400), + p = c(0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.05, 0.07, 0.1, 0.15)) +set.seed(2343) +# simulaciones +simulaciones <- parametros |> + mutate(muestra = map2(n, p, ~ simular_muestras(50000, .x, .y) |> select(rep, k))) |> + unnest(muestra) +# calcular_cobertura +calcular_cobertura <- function(simulaciones, construir_intervalo){ + # nombre de función + intervalo_nombre <- substitute(construir_intervalo) |> as.character() + cobertura_tbl <- simulaciones |> + mutate(intervalo = construir_intervalo(n, k)) |> + pull(intervalo) |> + bind_cols(simulaciones) |> + mutate(cubre = p >= inf & p <= sup) |> + group_by(n, p) |> + summarise(cobertura = mean(cubre), long_media = mean(sup - inf)) + cobertura_tbl |> mutate(tipo = intervalo_nombre) +} +``` + +```{r} +cobertura_wald <- calcular_cobertura(simulaciones, intervalo_wald) +cobertura_wald +``` + + +```{r} +graficar_cobertura <- function(cobertura_tbl){ + ggplot(cobertura_tbl, aes(x = p, y = cobertura, colour = tipo)) + + geom_hline(yintercept = 0.95, colour = "black") + + geom_line() + geom_point() + + facet_wrap(~n) + + ylim(0, 1) +} +cobertura_wald |> + graficar_cobertura() +``` + +La cobertura real es mucho más baja que la nominal en muchos casos, especialmente +cuando la $p$ es baja y $n$ es chica. Pero incluso para muestras relativamente grandes (100), +la cobertura es mala si $p$ es chica. + +Ahora probamos nuestro método alternativo: + +```{r} +cobertura_bayes <- calcular_cobertura(simulaciones, intervalo_bayes) +``` + + +```{r} +bind_rows(cobertura_wald, cobertura_bayes) |> + mutate(tipo = factor(tipo, levels = c('intervalo_wald', 'intervalo_bayes'))) |> + graficar_cobertura() +``` + +Y vemos que en general el intervalo de Bayes es superior al de Wald, en sentido +de que su cobertura real es más cercana a la nominal. El caso donde fallan los dos +es para muestras muy chicas $n=5, 10$, con probabilidades de éxito chicas $p\leq 0.02$. + +- Sin embargo, si tenemos información previa acerca del tamaño de la proporción que estamos +estimando, es posible obtener buena calibración con el método bayesiano. + +En este caso particular, **tenemos argumentos frecuentistas** +para utilizar el método bayesiano. Por ejemplo, si el INEGI utilizara estos +intervalos creíbles, un análisis de calibración de este tipo sostendría esa decisión. + +## Intervalos de Agresti-Coull {-} + +Un método intermedio que se usa para obtener mejores intervalos cuando +estimamos proporciones es el siguiente: + +- Agregar dos 1's y dos 0's a los datos. +- Utilizar el método de Wald con estos datos modificados. + +```{r} +intervalo_agresti_coull <- function(n, k){ + p_hat <- (k + 2)/ (n + 4) + ee_hat <- sqrt(p_hat * (1 - p_hat) / n) + tibble(inf = p_hat - 2 * ee_hat, sup = p_hat + 2 * ee_hat) +} +cobertura_ac <- calcular_cobertura(simulaciones, intervalo_agresti_coull) +``` + +```{r} +bind_rows(cobertura_wald, cobertura_bayes, cobertura_ac) |> + mutate(tipo = factor(tipo, levels = c('intervalo_wald', 'intervalo_bayes', 'intervalo_agresti_coull'))) |> + graficar_cobertura() +``` + +Que tiende a ser demasiado conservador para proporciones chicas: + +```{r} +graficar_cobertura(cobertura_ac) + + ylim(c(0.9, 1)) +``` + +**Conclusión 1**: Los intervalos de Agresti-Coull son una buena alternativa +para estimar proporciones como sustituto de los intervalos clásicos de Wald, aunque +tienden a ser muy conservadores para muestras chicas + +Idealmente podemos utilizar un método bayesiano +pues normalmente tenemos información inicial acerca de las proporciones que +queremos estimar. + +## Incorporando información inicial {-} + +Nótese que generalmente tenemos información acerca de la cantidad que +queremos estimar: por ejemplo, que proporción de visitantes de un sitio web +compra algo (usualmente muy baja, menos de 2%), qué proporción de personas tiene diabetes tipo 1 +(una proporción muy baja, menos de 1 por millar), o qué proporción de hogares tienen ingresos trimestrales +mayores a 150 mil pesos (menos de %5 con alta probabilidad). + +En este caso, tenemos que ajustar nuestra inicial. Por ejemplo, para el problema +de ingresos, podríamos usar una $\mathsf{Beta}(2, 100)$, cuyos cuantiles son: + +```{r} +# uno de cada 100 +a <- 2 +b <- 100 +beta_sims <- rbeta(5000, a, b) +quantile(beta_sims, c(0.01, 0.05, 0.50, 0.90, 0.99)) |> round(3) +``` +```{r} +qplot(beta_sims) +``` + +Veamos cómo se ven los intervalos bayesianos producidos con esta inicial: + +```{r} +crear_intervalo_bayes <- function(a, b){ + + intervalo_fun <- function(n, k){ + a_post <- k + a + b_post <- n - k + b + tibble(inf = qbeta(0.025, a_post, b_post), sup = qbeta(0.975, a_post, b_post)) + } + intervalo_fun +} +intervalo_bayes_2 <- crear_intervalo_bayes(a, b) +``` + + +```{r, cache = TRUE} +cobertura_bayes <- calcular_cobertura(simulaciones, + intervalo_bayes_2) +``` + + +```{r} +graficar_cobertura(bind_rows(cobertura_bayes, cobertura_ac) |> + filter(p < 0.05)) + + ylim(c(0.5, 1)) +``` + +Y vemos que la calibración es similar. Notemos sin embargo que +la longitud del del intervalo bayesiano es **mucho menor** que el de +Agresti-Coull cuando la muestra es chica: + +```{r} +ggplot(bind_rows(cobertura_bayes, cobertura_ac), + aes(x = p, y = long_media, colour = tipo)) + + geom_point() + facet_wrap(~n) +``` + +Cuando la muestra es chica, los intervalos de bayes son similares +a los iniciales, y mucho más cortos que los de Agresti-Coull. +Para muestras intermedias (50-100) los intervalos bayesianos +son más informativos que los de Agresti-Coull, con calibración similar, y +representan aprendizaje por encima de lo que sabíamos en la inicial. +Para muestras grandes, obtenemos resultados simililares. + +Por ejemplo: + +```{r} +set.seed(2131) +k <- rbinom(1, 50, 0.03) +k +intervalo_agresti_coull(50, k) |> round(3) +``` +es un intervalo muy grande que puede incluir valores negativos. En contraste, el intervalo +bayesiano es: + +```{r} +intervalo_bayes_2(50, k) |> round(3) +``` +Aún quitando valores negativos, los intervalos de Agresti-Coull son mucho más anchos. +La aproximación bayesiana, entonces, utiliza información previa +para dar un resultado considerablemente +más informativo, con calibración similar a Agresti-Coull. + +¿Aprendimos algo? Comparemos la posterior con la inicial: + +```{r} +beta_sims_inicial <- tibble(prop = rbeta(5000, a, b), dist = "inicial") +beta_sims_posterior <- tibble(prop = rbeta(5000, a + k, b + 50), dist = "posterior") +bind_rows(beta_sims_inicial, beta_sims_posterior) |> + ggplot(aes(x = prop, fill = dist)) + + geom_histogram(alpha = 0.5, position = "identity") +``` + +Donde vemos que no aprendimos mucho en este caso, pero nuestras creencias +sí cambiaron en comparación con la inicial. + +**Conclusión 2**: con el enfoque bayesiano podemos obtener intervalos +informativos con calibración razonable, incluso con información inicial que +no es muy precisa. Los intervalos de Agresti-Coull son poco informativos +para muestras chicas y/o proporciones chicas. + + +### Ejemplo: porporción de hogares de ingresos grandes {-} + +Usaremos los datos de ENIGH como ejemplo (ignorando el diseño, pero es posible hacer todas +las estimaciones correctamente) para estimar +el porcentaje de hogares que tienen ingreso corriente de más +de 150 mil pesos al trimestre. Suponemos que la muestra del enigh es la población, +y tomaremos una muestra iid de esta población. +Usamos la misma inicial que mostramos arriba, que es +una Beta con parámetros + +```{r} +c(a,b) +``` + +```{r} +set.seed(2521) +muestra_enigh <- read_csv("data/conjunto_de_datos_enigh_2018_ns_csv/conjunto_de_datos_concentradohogar_enigh_2018_ns/conjunto_de_datos/conjunto_de_datos_concentradohogar_enigh_2018_ns.csv") |> + select(ing_cor) |> + sample_n(120) |> + mutate(mas_150mil = ing_cor > 150000) +``` + +Un intervalo de 95% es entonces + +```{r} +k <- sum(muestra_enigh$mas_150mil) +k +intervalo_bayes_2(120, sum(muestra_enigh$mas_150mil)) |> round(3) +``` + +La media posterior es + +```{r} +prop_post <- (a + k) / (120 + b) +prop_post +``` + +¿Cuál es la verdadera proporción? + +```{r} +read_csv("data/conjunto_de_datos_enigh_2018_ns_csv/conjunto_de_datos_concentradohogar_enigh_2018_ns/conjunto_de_datos/conjunto_de_datos_concentradohogar_enigh_2018_ns.csv") |> + select(ing_cor) |> + mutate(mas_150mil = ing_cor > 150000) |> + summarise(prop_pob = mean(mas_150mil)) +``` +En este caso, nuestro intervalo cubre a la proporción poblacional. + +## Inferencia bayesiana y regularización {-} + +Como hemos visto en análisis y modelos anteriores, la posterior que usamos +para hacer inferencia combina aspectos de la inicial con la verosimilitud (los datos). +Una manera de ver esta combinación y sus beneficios es pensando en término de +*regularización* de estimaciones. + +- En las muestras hay variación. Algunas muestras particulares nos dan estimaciones +de máxima verosimilitud pobres de los parámetros de interés (estimaciones ruidosas). +- Cuando esas estimaciones pobres están en una zona de baja probabilidad de la +inicial, la estimación posterior tiende a moverse (o *encogerse*) hacia las +zonas de alta probabilidad de la inicial. +- Esto *filtra* ruido en las estimaciones. +- El mecanismo resulta en una reducción del **error cuadrático medio**, mediante +una reducción de la varianza de los estimadores (aunque quizá el sesgo aumente). + +Esta es una técnica poderosa, especialmente para problemas complejos donde tenemos +pocos datos para cada parámetro. En general, excluímos resultados que no concuerdan +con el conocimiento previo, y esto resulta en **mayor precisión** en las estimaciones. + +## Ejemplo: modelo normal y estaturas {-} + +Haremos un experimento donde simularemos muestras de los datos de cantantes. Usaremos el +modelo normal-gamma inverso que discutimos anteriormente, con la información inicial +que elicitamos. ¿Cómo se compara la estimación de máxima verosimilitud con la media +posterior? + +```{r} +# inicial para media, ver sección anterior para discusión (normal) +mu_0 <- 175 +n_0 <- 5 +# inicial para sigma^2 (gamma inversa) +a <- 3 +b <- 140 +``` + +Para este ejemplo chico, usaremos muestras de tamaño 5: + +```{r} +set.seed(3413) +# ver sección anterior para explicación de esta función +calcular_pars_posterior <- function(x, pars_inicial){ + # iniciales + mu_0 <- pars_inicial[1] + n_0 <- pars_inicial[2] + a_0 <- pars_inicial[3] + b_0 <- pars_inicial[4] + # muestra + n <- length(x) + media <- mean(x) + S2 <- sum((x - media)^2) + # sigma post + a_1 <- a_0 + 0.5 * n + b_1 <- b_0 + 0.5 * S2 + 0.5 * (n * n_0) / (n + n_0) * (media - mu_0)^2 + # posterior mu + mu_1 <- (n_0 * mu_0 + n * media) / (n + n_0) + n_1 <- n + n_0 + c(mu_1, n_1, a_1, b_1) +} +``` + +Y también de la sección anterior: + +```{r} +sim_params <- function(m, pars){ + mu_0 <- pars[1] + n_0 <- pars[2] + a <- pars[3] + b <- pars[4] + # simular sigmas + 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 +} +``` + + +```{r} +# simular muestras y calcular medias posteriores +simular_muestra <- function(rep, mu_0, n_0, a_0, b_0){ + cantantes <- lattice::singer |> + mutate(estatura_cm = 2.54 * height) |> + filter(str_detect(voice.part, "Tenor")) |> + sample_n(5, replace = FALSE) + pars_posterior <- calcular_pars_posterior(cantantes$estatura_cm, + c(mu_0, n_0, a_0, b_0)) + medias_post <- + sim_params(1000, pars_posterior) |> + summarise(across(everything(), mean)) |> + select(mu, sigma) + media <- mean(cantantes$estatura_cm) + est_mv <- c("mu" = media, + "sigma" = sqrt(mean((cantantes$estatura_cm - media)^2))) + bind_rows(medias_post, est_mv) |> + mutate(rep = rep, tipo = c("media_post", "max_verosim")) |> + pivot_longer(mu:sigma, names_to = "parametro", values_to = "estimador") +} +``` + + +```{r} +poblacion <- lattice::singer |> + mutate(estatura_cm = 2.54 * height) |> + filter(str_detect(voice.part, "Tenor")) |> + summarise(mu = mean(estatura_cm), sigma = sd(estatura_cm)) |> + pivot_longer(mu:sigma, names_to = "parametro", values_to = "valor_pob") +``` + +```{r, cache = TRUE} +errores <- map(1:2000, ~ simular_muestra(.x, mu_0, n_0, a, b)) |> + bind_rows() |> left_join(poblacion) |> + mutate(error = (estimador - valor_pob)) +ggplot(errores, aes(x = error, fill = tipo)) + + geom_histogram(bins = 20, position = "identity", alpha = 0.5) + facet_wrap(~parametro) +``` +Vemos claramente que la estimación de la desviación estándar de +nuestro modelo es claramente superior a la de máxima verosimilitud. En resumen: +```{r} +errores |> + group_by(tipo, parametro) |> + summarise(recm = sqrt(mean(error^2)) |> round(2)) |> + arrange(parametro) +``` + +Obtenemos una ganancia considerable en cuanto a la estimación de la desviación +estandar de esta población. Los estimadores de la media posterior +son superiores a los de máxima verosimilitud en términos de error cuadrático medio. + +Podemos graficar las dos estimaciones, muestra a muestra, para entender +cómo sucede esto: + +```{r} +errores |> + select(-error) |> + pivot_wider(names_from = tipo, values_from = estimador) |> + filter(parametro == "sigma") |> +ggplot(aes(x = max_verosim, y = media_post)) + + geom_abline(colour = "red") + + geom_hline(yintercept = sqrt(b/(a - 1)), lty = 2, color = 'black') + + geom_point() + + labs(subtitle = "Estimación de sigma") + + xlab("Estimador MV de sigma") + + ylab("Media posterior de sigma") + + coord_fixed() + + geom_segment(aes(x = 13, y = 11, xend = 13, yend = sqrt(b/(a - 1))), + colour='red', size=1, arrow =arrow(length = unit(0.5, "cm"))) + + geom_segment(aes(x = .5, y = 6, xend = .5, yend = sqrt(b/(a - 1))), + colour='red', size=1, arrow =arrow(length = unit(0.5, "cm"))) +``` + +Nótese como estimaciones demasiado bajas o demasiada altas son contraídas +hacia valores más consistentes con la inicial, lo cual resulta en +menor error. El valor esperado de $\sigma$ bajo la distribución inicial +se muestra como una horizontal punteada. + +## Ejemplo: estimación de proporciones {-} + +Ahora repetimos el ejercicio + +```{r} +# inicial +a <- 2 +b <- 100 +qbeta(c(0.01, 0.99), a, b) + +# datos +datos <- read_csv("data/conjunto_de_datos_enigh_2018_ns_csv/conjunto_de_datos_concentradohogar_enigh_2018_ns/conjunto_de_datos/conjunto_de_datos_concentradohogar_enigh_2018_ns.csv") |> + select(ing_cor) +# estimaciones +obtener_estimados <- function(datos){ + muestra_enigh <- datos |> + sample_n(120) |> + mutate(mas_150mil = ing_cor > 150000) + k <- sum(muestra_enigh$mas_150mil) + tibble(k = k, est_mv = k/120, media_post = (a + k) / (120 + b), pob = 0.02769) +} +estimadores_sim <- map(1:200, ~obtener_estimados(datos)) |> + bind_rows() +# calculo de errores +error_cm <- estimadores_sim |> + summarise(error_mv = sqrt(mean((est_mv - pob)^2)), + error_post = sqrt(mean((media_post - pob)^2))) +error_cm +``` +Podemos ver claramente que las medias posteriores están encogidas hacia valores +más chicos (donde la inicial tiene densidad alta) comparadas con las +estimaciones de máxima verosimilitud: + +```{r} +estimadores_sim_ag <- estimadores_sim |> + group_by(k, est_mv, media_post) |> + summarise(n = n()) +ggplot(estimadores_sim_ag, aes(x = est_mv, media_post, size = n)) + geom_point() + + geom_abline() +``` + +## Teoría de decisión {-} + +En esta parte (que sigue a @Wasserman a grandes rasgos), +discutimos brevemente teoría general que nos sirve para seleccionar +estimadores puntuales, y que esperemos ponga en contexto la parte +anterior que acabamos de discutir. Usaremos algunos conceptos que vimos +en la parte de propiedades de estimadores de máxima verosimilitud. + +Definimos una **función de pérdida** $L(\theta, \hat{\theta}_n)$, que +mide el costo de la discrepancia entre nuestro estimador + +$$\hat{\theta}_n = t(X_1,\ldots, X_n) = t(X)$$ + +y el verdadero valor $\theta$ +poblacional. Es posible considerar distintas funciones de pérdida, pero +como en secciones anteriores, usaremos la **pérdida cuadrática**, definida por: + +$$L(\theta, \hat{\theta}_n) = (\theta - \hat{\theta}_n)^2$$ + +Esta función toma distintos valores dependiendo de la muestra y del parámetro +$\theta$, y necesitamos resumirla para dar una evaluación de qué tan bueno +es el estimador $\hat{\theta}_n$. + +Ahora que hemos considerado tanto estadística bayesiana como frecuentista, +podemos pensar en resumir esta función de distintas maneras. + +Comenzamos pensando de manera frecuentista. En este caso, consideramos a $\theta$ como +un valor fijo, y nos preguntamos qué pasaría con la pérdida con distintas muestras +potenciales que podríamos observar. + +Definimos como antes el **riesgo** (frecuentista) como: + +$$R(\theta, t) = \mathbb{E}_X\left[ (\theta - \hat{\theta}_n)^2 \, \big| \, \theta\right]$$ + +donde promediamos a lo largo de las muestras posibles, con $\theta$ fijo. Esta +cantidad no nos dice necesariamente cómo escoger un buen estimador para +$\theta$, pues dependiendo de dónde está $\theta$ puede tomar valores distintos. + + + +Ahora vamos a pensar de manera bayesiana: en este caso, los datos serán fijos +una vez que los obervemos de manera que $\hat{\theta}_n$ está fijo, y el parámetro +$\theta$ es una cantidad aleatoria con distribución inicial $p(\theta)$. +Entonces consideraríamos el promedio sobre la posterior dado por: + + +$$\rho(t, X) = \mathbb{E}_{p(\theta|X)}\left[(\theta - \hat{\theta})^2 \, \big | \, X \right]$$ +que llamamos **riesgo posterior**. Esta cantidad se calcula con la posterior de +los parámetros dados los datos, y nos dice, una vez que vimos los datos, cómo es +el error de nuestro estimador. Nótese que esta cantidad no es útil para escoger +un estimador bueno $t$ antes de ver los datos, pero nos sirve para evaluar a un +estimador dados los datos. + +En el primer caso, promediamos sobre posibles muestras, y en el segundo por +valores posibles de $\theta$ para una muestra dada. + +### Ejemplo: riesgo frecuentista {-} + +Para observaciones bernoulli, el estimador de máxima verosimilitud es $\hat{p}_1 += k /n$, donde $n$ es el tamaño de muestra y $k$ el número de éxitos observados. +Podemos usar también como estimador la media posterior de un modelo +Beta-Bernoulli con inicial $a=2, b=2$, que nos daría $\hat{p}_2 = \frac{k + 2}{n ++ 4}$. + +Aunque podemos hacer los cálculos analíticos, aproximaremos el riesgo bajo +el error cuadrático usando simulación + +```{r} +perdida_cuad <- function(p, p_hat){ + (p - p_hat)^2 +} +# dos estimadores +t_1 <- function(n, x) x / n +t_2 <- function(n, x) (x + 2) / (n + 4) +estimar_riesgo <- function(n = 20, theta, perdida, reps = 10000){ + x <- rbinom(reps, n, theta) + # calcular estimaciones + theta_1 <- t_1(n, x) + theta_2 <- t_2(n, x) + # calcular pérdida + media_perdida <- tibble( + n = n, theta = theta, + estimador = c("MLE", "Posterior"), + riesgo = c(mean(perdida(theta, theta_1)), mean(perdida(theta, theta_2)))) + media_perdida +} +estimar_riesgo(n = 20, theta = 0.1, perdida = perdida_cuad) +``` + +Como dijimos, esta cantidad depende de $\theta$ que no conocemos. Así que +calculamos para cada valor de $\theta:$ + +Las funciones de riesgo $R(\theta, t_1)$ y $R(\theta, t_2)$ (dependen de $\theta$) +se ven aproximadamente como sigue: + +```{r} +p_seq <- seq(0, 1, 0.001) +riesgo_tbl <- map(p_seq, ~ estimar_riesgo(n = 20, theta = .x, perdida = perdida_cuad)) |> + bind_rows() +ggplot(riesgo_tbl, aes(x = theta, y = riesgo, colour = estimador)) + + geom_line() +``` + +Y vemos que el riesgo depende del verdadero valor del parametro: +en los extremos, el estimador de máxima verosimilitud tiene menos riesgo, +pero en el centro tiene más (esto es independiente del tipo de intervalos +que construyamos y su cobertura). + +La razón es que las estimaciones de tipo Agresti-Coull ($\theta_2$) +están contraídas hacia 0.5 (agregamos dos éxitos y dos fracasos). Esto produce +sesgo en la estimación para valores extremos de $\theta$. Sin embargo, +para valores centrales de $\theta$ tiene menos variabilidad +(por regularización) que el estimador de máxima verosimilitud, y sufre poco +de sesgo. + +### Ejemplo: riesgo posterior {-} + +Supongamos que la inicial es $\theta \sim \mathsf{Beta}(5,3)$ + +```{r} +estimar_riesgo_post <- function(n = 20, x, perdida, reps = 20000){ + # calcular estimaciones + theta_1 <- t_1(n, x) + theta_2 <- t_2(n, x) + # simular de posterior + theta_post <- rbeta(reps, x + 5, n - x + 3) + # calcular pérdida + media_perdida <- tibble( + n = n, x = x, + estimador = c("MLE", "Posterior"), + riesgo_post= c(mean(perdida(theta_post, theta_1)), mean(perdida(theta_post, theta_2)))) + media_perdida +} +estimar_riesgo_post(n = 20, x = 8, perdida = perdida_cuad) +``` +Como dijimos, esta cantidad depende de los datos $x$ que no hemos observado. Así que +calculamos para cada valor de $x$: + +Las funciones de pérdida promedio +$\rho(x, t_1)$ y $\rho(x, t_2)$ (dependen de $x$) +se ven aproximadamente como sigue: + +```{r} +x_seq <- seq(0, 20, 1) +riesgo_post_tbl <- + map(x_seq, ~ estimar_riesgo_post(n = 20, x = .x, perdida = perdida_cuad)) |> + bind_rows() +ggplot(riesgo_post_tbl, aes(x = x, y = riesgo_post, colour = estimador)) + + geom_line() + geom_point() +``` + +Donde vemos que la pérdida del estimador bayesiano es mejor para valores extremos +de número de éxitos observado $x$, pero tiene más riesgo posterior +para valores chicos de $x$. En general es mejor el estimador $\theta_2$. El +estimador de máxima verosimilitud tiene más riesgo en los extremos, lo que esperaríamos +porque no tenemos la regularización que aporta la posterior. Igualmente, +vemos más riesgo para valores chicos de $x$ que para valores grandes: esto +es porque la inicial está concentrada en valores reslativamente grandes +de $\theta$. + + + +## Riesgo de Bayes {-} + +Finalmente, podemos crear un resumen unificado considerando: + +- Si no conocemos el valor del parámetro $\theta$, podemos promediar +el riesgo frecuentista con la inicial $p(\theta)$ +- Si no conocemos los datos observados, podemos promediar usando +datos generados por la marginal $p(x)$ de $x$ bajo el modelo de datos +$p(x|\theta)$ y la inicial $p(\theta)$. + +Por la ley de la esperanza iterada, estos dos resultados son iguales. La cantidad +resultante + +$$r(t) = \int R(\theta,t) p(\theta)\, d\theta = \int r(x, t)p(x|\theta)p(\theta)\, d\theta\, dx$$ +Se llama **riesgo de Bayes** para el estimador $t$. + +### Ejemplo {-} + +Podemos calcular + +```{r} +marginal_tbl <- function(n = 20, m = 5000){ + theta <- rbeta(m, 5, 3) + x <- rbinom(m, size = n, p = theta) + tibble(x = x) |> group_by(x) |> + summarise(n_x = n()) +} +riesgo_post_tbl |> left_join(marginal_tbl()) |> + group_by(estimador) |> + summarise(riesgo_bayes = sum(riesgo_post * n_x) / sum(n_x)) +``` + +o también + +```{r} +theta_tbl <- tibble(theta = rbeta(50000, 5, 3) |> round(3)) |> + group_by(theta) |> + summarise(n_x = n()) + +riesgo_tbl |> left_join(theta_tbl) |> + mutate(n_x = ifelse(is.na(n_x), 0, n_x)) |> + group_by(estimador) |> + summarise(riesgo_bayes = sum(riesgo * n_x) / sum(n_x)) +``` + + + + +Ahora consideremos cómo decidiríamos, desde el punto de vista Bayesiano, qué +estimador usar: + +1. (Estimador de Bayes) Si tenemos los datos $X$, escogeríamos una +función $t_X$ que minimice el riesgo posterior $\rho(t, X)$, y nuestro +estimador es $\hat{\theta}_n = t_X (X)$. +2. (Regla de Bayes) Si no tenemos los datos, escogeríamos el estimador una función +$t$ que minimice el riesgo de Bayes $r(t)$, y estimaríamos usando $\hat{\theta}_n = t(X)$ +3. Pero como el riesgo de Bayes es el promedio del riesgo posterior, la +solución de 1 nos da la solución de 2. Es decir, el estimador que escogemos +condicional a los datos $X$ es el mismo que escogeríamos antes de escoger los +datos, dada una distribución inicial $p(\theta).$ + + +Por ejemplo, es posible demostrar que bajo la pérdida cuadrática, la regla de +Bayes es utilizar la media posterior, bajo la pérdida absoluta, la mediana +posterior, etc. + +Este estimador de Bayes tiene sentido desde el punto de vista frecuentista, +también, porque **minimiza el riesgo frecuentista promedio**, suponiendo la +inicial $p(\theta)$. Por ejemplo, para la pérdida cuadrática podemos usar la +descomposición de sesgo y varianza y obtenemos: + +$$r(t) = \mathbb{E}[R(\theta,t)] = \mathbb{E}[ \mathsf{Sesgo}_\theta^2(t)] +\mathbb{E}[\mathsf{Var}_\theta(t)] $$ + +Podemos ver entonces que el estimador de Bayes, en este caso la media posterior, +es resultado de minimizar la suma de estas dos cantidades: por eso puede +incurrir en sesgo, si ese costo se subsana con una reducción considerable de la +varianza. Los estimadores insesgados que vimos en esta sección fueron subóptimos +en muchos casos justamente porque son insesgados, e incurren en varianza grande. + +```{block2, type = 'ejercicio'} +Regresa a los ejemplos anteriores donde evaluamos el desempeño de la media +posterior en varios ejemplos. Muestra en las gráficas dónde ves el balance +entre sesgo y varianza que cumplen cuando los comparamos con estimadores +insesgados. +``` + + +Desde el punto de vista frecuentista, la cuestión es más complicada y +hay varias maneras de proceder. En primer lugar, comenzaríamos con el +riesgo frecuentista $R(\theta, t)$. Una idea es,por ejemplo, calcular el +riesgo máximo: +$$R_{\max} (t) = \underset{\theta}{\max} R(\theta, t).$$ +En nuestro ejemplo de arriba el máximo se alcanza en 0.5, y tomaríamos +eso evaluación de los estimadores $\theta_1$ o $\theta_2$. Buscaríamos entonces +estimadores que minimicen este máximo, es decir, **estimadores minimax**. Pero también +es posible enfocar este problema considerando sólo estimadores *insesgados*, lo que +nos lleva por ejemplo a buscar estimadores con mínima varianza. También podemos +enfocarnos en buscar estimador *admisibles*, que son aquellos cuyo riesgo no está +dominado para toda $\theta$ por otro estimador, y así sucesivamente. + +Finalmente, es posible demostrar (ver @Wasserman) que típicamente, para muestras +grandes, el estimador de máxima verosimilitud es cercano +a ser minimax y además es una regla de Bayes. Estas son buenas propiedades, pero +debemos contar con que el régimen asintótico se cumpla aproximadamente. + + + diff --git a/16-bayes-mcmc.Rmd b/16-bayes-mcmc.Rmd new file mode 100644 index 0000000..92a75e2 --- /dev/null +++ b/16-bayes-mcmc.Rmd @@ -0,0 +1,1771 @@ +# Métodos de Cadenas de Markov Monte Carlo + + +```{r setup, include=FALSE, message=FALSE} +library(tidyverse) +library(patchwork) +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center', fig.width = 5, fig.height=3, cache = TRUE) +comma <- function(x) format(x, digits = 2, big.mark = ",") +theme_set(theme_minimal()) +``` + +Hasta ahora, hemos considerado modelos bayesianos *conjugados*, donde la +posterior tiene una forma conocida. Esto nos permitió simular directamente +de la posterior usando las rutinas estándar de `R`, o utilizar cálculos teóricos +o funciones estándar de `R` para calcular resúmenes de interés, como medias o +medianas posteriores o intervalos de credibilidad. + +Sin embargo, en aplicaciones rara vez es factible este tipo de análisis tan +simple, pues: + +1. Los modelos que estamos considerando son más complejos y la distribución posterior +conjunta de los parámetros no tiene una forma simple conocida. +2. Queremos usar distribuciones iniciales que no son conjugadas para utilizar correctamente +nuestra información inicial. + +Recordamos que tenemos expresiones explícitas para la inicial $p(\theta)$ y la verosimilitud +$p(x|\theta)$, así que conocemos explícitamente la posterior, módulo la constante de +normalización, + +$$p(\theta|x) \propto p(x|\theta) \, p(\theta).$$ + + +Supongamos por ejemplo que quisiéramos calcular las medias posteriores de los +parámetros $\theta$. En teoría, tendríamos que calcular + +$$\hat \theta = \mathbb{E}[{\theta}\, |\, x] = \int \theta \, p(\theta|x) \, d\theta$$ +Entonces es necesario calcular también $p(x)$, que resulta de la integral +$$p(x) = \int p(x|\theta) \, p(\theta)\, d\theta$$ + +Si no tenemos expresiones analíticas simples, tendremos que aproximar numéricamente +estas integrales de alguna forma. + +1. Si la posterior tiene una forma conocida, podemos calcular cantidades de interés usando +fórmulas o rutinas de simulación de distribuciones conocidas que producen muestras independientes. + +Cuando la posterior no tiene una forma conocida, sin embargo: + +2. Podemos intentar usar integración numérica usual. Como veremos, este enfoque no es muy escalable. +3. Podemos usar simulaciones bajo cadenas de Markov (**Markov Chain Monte Carlo**, MCMC), +que es un enfoque más escalable. + +Mucho del uso generalizado actual de la estadística bayesiana se debe a que gracias al +poder de cómputo disponible y los métodos MCMC, no estamos restringidos al uso de 1 y 2, +que tienen desventajas grandes. Primero mostraremos cómo el método de integración +por subdivisión no es escalable. + + +## Integrales mediante subdivisiones {-} + +Como tenemos una expresión analítica para el integrando, podemos intentar una +rutina numérica de integración. Una vez calculada, podríamos entonces +usar otra rutina numérica para calcular las medias posteriores $\hat{\theta}$. + +Las rutinas usuales de integración pueden sernos útiles cuando el número de parámetros +es chico. Consideremos primero el caso de 1 dimensión, y supongamos que $a\leq\theta\leq b$. + +Si dividimos el rango de $\theta$ en intervalos determinados +por $a = \theta^1<\theta^2<\cdots \theta^M =b$, tales que $\Delta\theta = \theta^{i+1} -\theta^{i}$, +podríamos aproximar con + +$$p(x) \approx \sum_{i=1}^M p(x|\theta^i)p(\theta^i) \Delta\theta$$ +Lo que requiere $M$ evaluaciones del factor $p(x|\theta)p(\theta)$. Podríamos usar +por ejemplo $M=100$ para tener precisión razonable. + +### Ejemplo: estimación de una proporción {-} + +Teníamos que $p(S_n = k|\theta) \propto \theta^k(1-\theta)^{n-k}$ cuando observamos $k$ éxitos +en $n$ pruebas independientes. Supongamos que nuestra inicial es $p(\theta) = 2\theta$ +(checa que es una densidad), es decir, creemos que es más probable a priori +observar proporciones altas. Podemos integrar numéricamente + +```{r} +crear_log_post <- function(n, k){ + function(theta){ + verosim <- k * log(theta) + (n - k) * log(1 - theta) + inicial <- log(theta) + log_p_factor <- verosim + inicial + log_p_factor + } +} +# observamos 3 éxitos en 4 pruebas: +log_post <- crear_log_post(4, 3) +prob_post <- function(x) { exp(log_post(x))} +# integramos numéricamente +p_x <- integrate(prob_post, lower = 0, upper = 1, subdivisions = 100L) +p_x +``` + +Y ahora podemos calcular la media posterior: + +```{r} +media_funcion <- function(theta){ + theta * prob_post(theta) / p_x$value +} +integral_media <- integrate(media_funcion, lower = 0, upper = 1, subdivisions = 100L) +media_post <- integral_media$value +media_post +``` +Podemos verificar nuestro trabajo pues sabemos que la posterior es $\mathsf{Beta}(5, 2)$ +cuya media es +```{r} +5/(2+5) +``` + +Y podríamos intentar una estrategia similar, por ejemplo, para calcular intervalos +de credibilidad. Sin embargo, veremos abajo que este método no escala con el número de +parámetros. + +### Más de un parámetro {-} + +Ahora supongamos que tenemos $2$ parámetros. Dividiríamos cada parámetro +en 100 intervalos, y luego tendríamos que calcular + +$$p(x) \approx \sum_{i=1}^M \sum_{j=1}^M p(x|\theta_1^i, \theta_2^j)p(\theta_1^i, \theta_2^j) \Delta\theta_1\Delta\theta_2$$ +Y esto requeriría $M^2 = 10,000$ evaluaciones de $p(x|\theta)p(\theta)$. + +Si tenemos $p$ parámetros, entonces tendríamos que hacer $M^p$ evaluaciones de la +posterior. Incluso cuando $p=10$, **esta estrategia es infactible**, pues tendríamos +que hacer más de millones de millones de millones de evaluaciones de la posterior. Si sólo +tenemos esta técnica disponible, el análisis bayesiano está considerablemente +restringido. Regresión bayesiana con unas 10 covariables por ejemplo, no podría hacerse. + +De modo que tenemos que replantearnos cómo atacar el problema de calcular o aproximar +estas integrales. + +## Métodos Monte Carlo {-} + +En varias ocasiones anteriormente hemos usado el método Monte Carlo para +aproximar integrales: por ejemplo, para calcular medias posteriores. + +Supongamos que tenemos una densidad $p(\theta)$. + +```{block2, type='mathblock'} +**Integración Monte Carlo**. Supongamos que queremos calcular el valor esperado de +$g(X)$, donde $X\sim p(X\,|\,\theta).$ Es decir, la variable aleatoria $X$ se +distribuye de acuerdo al modelo probabilistico $p(X \, | \, \theta),$ de tal forma que +lo que nos interesa calcular es + +$$\mathbb{E}[g(X)] = \int g(x) p(x|\theta)\, dx.$$ + +Si tomamos una muestra +$x^{(1)},x^{(2)}, \ldots x^{(N)} \overset{iid}{\sim} p(x|\theta)$, entonces + +$$\mathbb{E}[g(X)] \approx \, \frac1N \, \sum_{n = 1}^N g(x^{(n)})$$ + +cuando $N$ es grande. +``` + +Esto es simplemente una manera de escribir la ley de los grandes números, y hemos +aplicado este teorema en varias ocasiones. Nos ha sido útil cuando +**sabemos cómo simular de distribución** $p(\theta | x)$ (usando alguna rutina de `R`, por +ejemplo, o usando un método estándar como inversión de la función de distribución acumulada). + +### Ejemplo {-} +En este ejemplo repetimos cosas que ya hemos visto. En el caso de estimación +de una proporción $\theta$, tenemos como inicial +$p(\theta) \propto \theta$, que es $\mathsf{Beta}(2,1)$. Si observamos 3 éxitos en 4 pruebas, +entonces sabemos que la posterior es $p(\theta|x)\propto \theta^4(1-\theta)$, que +es $\mathsf{Beta}(5, 2)$. Si queremos calcular media y segundo momento posterior para $\theta$, +en teoría necesitamos calcular + +$$\mu = \int_0^1 \theta p(\theta|X = 3)\, d\theta,\,\, \mu_2=\int_0^1 \theta^2 p(\theta|X = 3)\, d\theta$$ + +integramos con Monte Carlo + +```{r} +theta <- rbeta(10000, 5, 2) +media_post <- mean(theta) +momento_2_post <- mean(theta^2) +c(media_post, momento_2_post) +``` + +Y podemos aproximar de esta manera cualquier cantidad de interés que esté basada +en integrales, como probabilidades asociadas a $\theta$ o cuantiles asociados. +Por ejemplo, podemos aproximar fácilmente $P(e^{\theta}> 2|x)$ haciendo + +```{r} +mean(exp(theta) > 2) +``` +y así sucesivamente. + +Este enfoque, sin embargo, es mucho más flexible y poderoso. + +### Ejemplo: varias pruebas independientes {-} + +Supongamos que probamos el nivel de gusto para 4 sabores distintos de una paleta. Usamos +4 muestras de aproximadamente +50 personas diferentes para cada sabor, y cada uno evalúa si le gustó mucho o no. +Obtenemos los siguientes resultados: + +```{r} +datos <- tibble( + sabor = c("fresa", "limón", "mango", "guanábana"), + n = c(50, 45, 51, 50), gusto = c(36, 35, 42, 29)) %>% + mutate(prop_gust = gusto / n) +datos +``` + +Usaremos como inicial $\mathsf{Beta}(2, 1)$ (pues hemos obervado cierto sesgo de +cortesía en la calificación de sabores, y no es tan probable tener valores muy +bajos) para todos los sabores, es decir $p(\theta_i)$ es la funcion de densidad +de una $\mathsf{Beta}(2, 1)$. La inicial conjunta la definimos entonces, usando +idependiencia inicial, como + +$$p(\theta_1,\theta_2, \theta_3,\theta_4) = p(\theta_1)p(\theta_2)p(\theta_3)p(\theta_4).$$ +Pues inicialmente establecemos que ningún parámetro da información sobre otro: +saber que mango es muy gustado no nos dice nada acerca del gusto por fresa. Bajo +este supuesto, y el supuesto adicional de que las muestras de cada sabor son +independientes, podemos mostrar que las posteriores son independientes: + +$$p(\theta_1,\theta_2,\theta_3, \theta_4|k_1,k_2,k_3,k_4) = p(\theta_4|k_1)p(\theta_4|k_2)p(\theta_4|k_3)p(\theta_4|k_4)$$ + +De forma que podemos trabajar individualmente con cada muestra. Calculamos los parámetros de las posteriores individuales: + +```{r} +datos <- datos %>% + mutate(a_post = gusto + 2, b_post = n - gusto + 1) +datos +``` + +Ahora nos preguntamos, ¿cuál es la probabilidad posterior de que mango sea el sabor +más preferido de la población? Conocemos la posterior para cada parámetro, y sabemos +que los parámetros son independientes para la posterior. Eso quiere decir +que podemos simular de cada parámetro independientemente para obtener simulaciones +de la conjunta posterior. + +```{r} +simular_conjunta <- function(rep, datos){ + datos %>% mutate(valor_sim = map2_dbl(a_post, b_post, ~ rbeta(1, .x, .y))) %>% + select(sabor, valor_sim) +} +simular_conjunta(1, datos) +``` + +```{r} +# esta no es una manera muy rápida, podríamos calcular todas las +# simulaciones de cada parámetro de manera vectorizada +sims_posterior <- tibble(rep = 1:5000) %>% + mutate(sims = map(rep, ~ simular_conjunta(.x, datos))) %>% + unnest(cols = sims) +sims_posterior +``` + + +Y ahora podemos aproximar fácilmente la probabilidad de interés: + +```{r} +sims_posterior %>% + group_by(rep) %>% + mutate(sabor = sabor[which.max(valor_sim)]) %>% + group_by(sabor) %>% + count() %>% + ungroup() %>% + mutate(prop = n / sum(n)) +``` +Y vemos que los mejores sabores son mango y limón. La probabilidad posterior de +que mango sea el sabor preferido por la población es de 66%. La integral correspondiente +no es trivial. + + +```{block2, type='ejercicio'} +- ¿Cuáles son las probabilidades a priori de que cada sabor sea el preferido +por la población? +- ¿Cuál es la integral correspondiente a las probabilidades que acabamos de calcular? + ¿Qué tan fácil es hacer esta integral de manera analítica? +- Calcula la probabilidad de que mango sea preferido a limón? +- ¿Qué conclusión práctica sacas de estos resultados? +``` + + +## Simulando de la posterior {-} + +Hemos establecido que podemos contestar varias preguntas de inferencia +usando simulación Monte Carlo, y que este enfoque es potencialmente +escalable (en contraste con métodos de integración numérica por cuadrícula). Ahora +el problema que necesitamos resolver es el siguiente: + +- Conocemos $p(\theta |x)$ módulo una constante de integración. +- En general, $p(\theta|x)$ no tiene una forma reconocible que corresponda a un +simulador estándar. +- ¿Cómo simulamos de esta posterior cuando sólo sabemos calcular $p(x|\theta)p(\theta)$? + +Hay varias maneras de hacer esto. Presentaremos los algoritmos en términos +de una distribución cualquiera $p(\theta) = K f(\theta)$, donde sólo conocemos +la función $f(\theta)$. + + +## Ejemplo de islas + +Comenzamos revisando el ejemplo de las islas en @Kruschke (7.2) para tener más intuición de cómo funciona este algoritmo. + +## ¿Por qué funciona Metrópolis? {-} + +Supongamos que un vendedor de *Yakult* trabaja a lo largo de una cadena de +islas: + +* Constantemente viaja entre las islas ofreciendo sus productos; + +* Al final de un día de trabajo decide si permanece en la misma isla o se +transporta a una de las $2$ islas vecinas; + +* El vendedor ignora la distribución de la población en las islas y el número +total de islas; sin embargo, una vez que se encuentra en una isla puede +investigar la población de la misma y también de la isla a la que se propone +viajar después. + +* El objetivo del vendedor es visitar las islas de manera proporcional a la +población de cada una. Con esto en mente el vendedor utiliza el siguiente +proceso: + 1) Lanza un volado, si el resultado es águila se propone ir a la isla +del lado izquierdo de su ubicación actual y si es sol a la del lado derecho. + 2) Si la isla propuesta en el paso anterior tiene población mayor a la +población de la isla actual, el vendedor decide viajar a ella. Si la isla vecina +tiene población menor, entonces visita la isla propuesta con una probabilidad que +depende de la población de las islas. Sea $P^*$ la población de la isla +propuesta y $P_{t}$ la población de la isla actual. Entonces el vendedor +cambia de isla con probabilidad +$$q_{mover}=P^*/P_{t}$$ + +A la larga, si el vendedor sigue la heurística anterior la probabilidad de que +el vendedor este en alguna de las islas coincide con la población relativa de +la isla. + +```{r, fig.height=6, fig.width=3.5} +islas <- tibble(islas = 1:10, pob = 1:10) +camina_isla <- function(i){ # i: isla actual + u <- runif(1) # volado + v <- ifelse(u < 0.5, i - 1, i + 1) # isla vecina (índice) + if (v < 1 | v > 10) { # si estás en los extremos y el volado indica salir + return(i) + } + p_move = ifelse(islas$pob[v] > islas$pob[i], 1, islas$pob[v] / islas$pob[i]) + u2 <- runif(1) + if (p_move > u2) { + return(v) # isla destino + } + else { + return(i) # me quedo en la misma isla + } +} +pasos <- 100000 +iteraciones <- numeric(pasos) +iteraciones[1] <- sample(1:10, 1) # isla inicial +for (j in 2:pasos) { + iteraciones[j] <- camina_isla(iteraciones[j - 1]) +} +caminata <- tibble(pasos = 1:pasos, isla = iteraciones) +plot_caminata <- ggplot(caminata[1:1000, ], aes(x = pasos, y = isla)) + + geom_point(size = 0.8) + + geom_path(alpha = 0.5) + + coord_flip() + + labs(title = "Caminata aleatoria") + + scale_y_continuous(expression(theta), breaks = 1:10) + + scale_x_continuous("Tiempo") +plot_dist <- ggplot(caminata, aes(x = isla)) + + geom_histogram() + + scale_x_continuous(expression(theta), breaks = 1:10) + + labs(title = "Distribución objetivo", + y = expression(P(theta))) +plot_caminata / plot_dist +``` + +Entonces: + +* Para aproximar la distribución objetivo debemos permitir que el vendedor +recorra las islas durante una sucesión larga de pasos y registramos sus visitas. + +* Nuestra aproximación de la distribución es justamente el registro de sus +visitas. + +* Más aún, debemos tener cuidado y excluir la porción de las visitas que se +encuentran bajo la influencia de la posición inicial. Esto es, debemos excluir +el **periodo de calentamiento**. + +* Una vez que tenemos un registro _largo_ de los viajes del vendedor (excluyendo +el calentamiento) podemos aproximar la distribución objetivo +simplemente contando el número relativo de veces que el vendedor visitó +dicha isla. + +```{r, warning=FALSE, message=FALSE, fig.width=8, fig.height=7.5} +t <- c(1:10, 20, 50, 100, 200, 1000, 5000) +plots_list <- map(t, function(i){ + ggplot(caminata[1:i, ], aes(x = isla)) + + geom_histogram() + + labs(y = "", x = "", title = paste("t = ", i, sep = "")) + + scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11)) +}) +wrap_plots(plots_list) +``` + + +Escribamos el algoritmo, para esto indexamos las islas por el valor +$\theta$, es así que la isla del extremo oeste corresponde a $\theta=1$ y la +población relativa de cada isla es $P(\theta)$: + +1. El vendedor se ubica en $\theta^{(i)}$ y propone moverse a la izquierda +o derecha con probabilidad $0.5$. +El rango de los posibles valores para moverse, y la probabilidad de proponer +cada uno se conoce como **distribución propuesta**, en nuestro ejemplo sólo +toma dos valores cada uno con probabilidad $0.5$. + +2. Una vez que se propone un movimiento, decidimos si aceptarlo. La decisión de +aceptar se basa en el valor de la distribución **objetivo** en la posición +propuesta, relativo al valor de la distribución objetivo en la posición actual: +$$\alpha=\min\bigg\{\frac{P(\theta^*)}{P(\theta^{(i)})},1\bigg\},$$ +donde $\alpha$ denota la probabilidad de hacer el cambio de isla. + +Notemos que la distribución objetivo $P(\theta)$ no necesita estar normalizada, +esto es porque lo que nos interesa es el cociente $P(\theta^*)/P(\theta^{(i)})$. + +3. Una vez que propusimos un movimiento y calculamos la probabilidad de aceptar +el movimiento aceptamos o rechazamos el movimiento generando un valor de una +distribución uniforme, si dicho valor es menor a la probabilidad de cambio, +$\alpha,$ entonces hacemos el movimiento. + +Entonces, para utilizar el algoritmo necesitamos ser capaces de: + +* Generar un valor de la distribución propuesta, que hemos denotado por $q,$ +(para crear $\theta^*$). + +* Evaluar la distribución objetivo en cualquier valor propuesto (para calcular +$P(\theta^*)/P(\theta^{(i)})$). + +* Generar un valor uniforme (para movernos con probabilidad $\alpha$). + +Las $3$ puntos anteriores nos permiten generar muestras aleatorias de la +distribución objetivo, sin importar si esta está normalizada. Esta técnica es +particularmente útil cuando cuando la distribución objetivo es una posterior +proporcional a $p(x|\theta)p(\theta)$. + + +Para entender porque funciona el algoritmo de Metrópolis hace falta entender $2$ +puntos, primero que la distribución objetivo es **estable**: si la probabilidad +_actual_ de ubicarse en una posición coincide con la probabilidad en la +distribución objetivo, entonces el algoritmo preserva las probabilidades. + +```{r, warning=FALSE, message=FALSE, fig.width=8, fig.height=7.5} +library(expm) +transMat <- function(P){ # recibe vector de probabilidades (o población) + T <- matrix(0, 10, 10) + n <- length(P - 1) # número de estados + for (j in 2:n - 1) { # llenamos por fila + T[j, j - 1] <- 0.5 * min(P[j - 1] / P[j], 1) + T[j, j] <- 0.5 * (1 - min(P[j - 1] / P[j], 1)) + + 0.5 * (1 - min(P[j + 1] / P[j], 1)) + T[j, j + 1] <- 0.5 * min(P[j + 1] / P[j], 1) + } + # faltan los casos j = 1 y j = n + T[1, 1] <- 0.5 + 0.5 * (1 - min(P[2] / P[1], 1)) + T[1, 2] <- 0.5 * min(P[2] / P[1], 1) + T[n, n] <- 0.5 + 0.5 * (1 - min(P[n - 1] / P[n], 1)) + T[n, n - 1] <- 0.5 * min(P[n - 1] / P[n], 1) + T +} +T <- transMat(islas$pob) +w <- c(0, 1, rep(0, 8)) +t <- c(1:10, 20, 50, 100, 200, 1000, 5000) +expT <- map_df(t, ~data.frame(t = ., w %*% (T %^% .))) +expT_long <- expT %>% + gather(theta, P, -t) %>% + mutate(theta = parse_number(theta)) +ggplot(expT_long, aes(x = theta, y = P)) + + geom_bar(stat = "identity", fill = "darkgray") + + facet_wrap(~ t) + + scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11)) +``` + +El segundo punto es que el proceso converge a la distribución objetivo. +Podemos ver, (en nuestro ejemplo sencillo) que sin importar el punto de inicio +se alcanza la distribución objetivo. + +```{r, warning=FALSE, message=FALSE, fig.width=8, fig.height=7.5} +inicio_p <- function(i){ + w <- rep(0, 10) + w[i] <- 1 + t <- c(1, 10, 50, 100) + exp_t <- map_df(t, ~ data.frame(t = .x, inicio = i, w %*% (T %^% .))) %>% + gather(theta, P, -t, -inicio) %>% + mutate(theta = parse_number(theta)) + exp_t +} +exp_t <- map_df(c(1, 3, 5, 9), inicio_p) +ggplot(exp_t, aes(x = as.numeric(theta), y = P)) + + geom_bar(stat = "identity", fill = "darkgray") + + facet_grid(inicio ~ t) + + scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11)) +``` + + + +## Método de Metrópolis {-} + +En el método de Metrópolis, uno de los más antiguos, comenzamos +con un valor inicial de los parámetros $\theta^{(0)}$ en el soporte de $p(\theta)$, +es decir $p(\theta^{(0)})>0.$ + +Para $i=1, \ldots, M$, hacemos: + +1. Partiendo de $\theta^{(i)}$, hacemos un salto +corto en una dirección al azar para obtener una propuesta $\theta^* \sim q(\theta \, |\, \theta^{(i)}).$ +2. Aceptamos or rechazamos el salto: + - Si $\alpha = \frac{f(\theta^*)}{f(\theta^{(i)})} \geq 1$, aceptamos el salto + y ponemos $\theta^{(i+1)}=\theta^*$. Regresamos a 1 para la siguiente + iteración $i\leftarrow i + 1.$ + - Si $\alpha = \frac{f(\theta^*)}{f(\theta^{(i)})} < 1$, entonces aceptamos + con probabilidad $\alpha$ el salto, ponemos $\theta^{(i+1)}=\theta^*$ y + regresamos a 1 para la siguiente iteración $i\leftarrow i + 1$. Si rechazamos + el salto, ponemos entonces $\theta^{(i+1)}=\theta^{(i)}$ y regresamos a 1 para + la siguiente iteración $i\leftarrow i + 1.$ + +Requerimos también que la función que propone los saltos sea simétrica: es +decir, $q(\theta^*|\theta^{(i)})$ debe ser igual a $q(\theta^{(i)}|\theta^*)$. +Se puede modificar el algoritmo para tratar con una propuesta que no sea +simétrica. + +Una elección común es escoger $q(\theta^* |\theta^{(i)})$ como +una $\mathsf{N}(\theta^{(i)}, \sigma_{salto})$. + +```{block,type="comentario"} +En este curso, escribiremos varios métodos de cadenas de Markov para estimación +Monte Carlo (*Markov Chain Monte Carlo*, MCMC) desde cero para entender los +básicos de cómo funciona. Sin embargo, **en la práctica no hacemos esto**, sino +que usamos software estándar (Stan, JAGS, BUGS, etc.) para hacer este trabajo. + +**Expertos** en MCMC, métodos numéricos, y estadística a veces escriben partes +de sus rutinas de simulación, y pueden lograr mejoras de desempeño +considerables. Excepto para modelos simples, esto no es trivial de hacer +garantizando resultados correctos. +``` + +En resumen, todo el código de esta sección es de carácter ilustrativo. +Utiliza implementaciones establecidas en las aplicaciones. + +Abajo implementamos el algoritmo con un salto de tipo normal: + +```{r} +crear_metropolis <- function(fun_log, sigma_salto = 0.1){ + # la entrada es la log posterior + iterar_metropolis <- function(theta_inicial, n){ + p <- length(theta_inicial) + nombres <- names(theta_inicial) + iteraciones <- matrix(0, nrow = n, ncol = p) + colnames(iteraciones) <- nombres + iteraciones[1,] <- theta_inicial + for(i in 2:n){ + theta <- iteraciones[i - 1, ] + theta_prop <- theta + rnorm(p, 0, sigma_salto) + # exp(log(p) - log(q)) = p/q + cociente <- exp(fun_log(theta_prop) - fun_log(theta)) + if(cociente >= 1 || runif(1,0,1) < cociente){ + iteraciones[i, ] <- theta_prop + } else { + iteraciones[i, ] <- theta + } + } + iteraciones_tbl <- iteraciones %>% + as_tibble() %>% + mutate(iter_num = row_number()) %>% + select(iter_num, everything()) + iteraciones_tbl + } + iterar_metropolis +} + +``` + +E intentamos simular de una exponencial no normalizada: + +```{r} +exp_no_norm <- function(x) { + z <- ifelse(x > 0, exp(-0.5 * x), 0) + log(z) +} + +iterador_metro <- crear_metropolis(exp_no_norm, sigma_salto = 0.25) +sims_tbl <- iterador_metro(c(theta = 0.5), 50000) +ggplot(sims_tbl, aes(x = theta)) + geom_histogram() +``` + +Ahora probemos con una $\mathsf{Beta}(3, 2):$ + +```{r} +beta_no_norm <- function(x) { + z <- ifelse(x > 0 && x < 1, (x^2)*(1-x), 0) + log(z) +} + +iterador_metro <- crear_metropolis(beta_no_norm, sigma_salto = 0.04) +sims_metro_tbl <- iterador_metro(c(theta = 0.5), 50000) +sims_indep_tbl <- tibble(iter_num = 1:30000, theta = rbeta(30000, 3, 2)) +g_1 <- ggplot(sims_metro_tbl, aes(x = theta)) + geom_histogram() + + labs(subtitle = "Metrópolis") +g_2 <- ggplot(sims_indep_tbl, aes(x = theta)) + + geom_histogram() + + labs(subtitle = "rbeta") +g_1 + g_2 +``` +Y vemos que esto funciona. + +Nótese sin embargo un aspecto de estas simulaciones que no habíamos encontrado +en el curso. Aunque la distribución final de las simulaciones es muy cercana +a la de la distribución que queremos simular, lo cual era nuestro propósito, +las **simulaciones no son extracciones independientes** de esa distribución. + +La construcción del algoritmo muestra eso, pero podemos también graficar las +simulaciones: + +```{r, fig.width = 8} +g_metropolis <- sims_metro_tbl %>% + filter(iter_num < 500) %>% + ggplot(aes(x = iter_num, y = theta)) + + geom_line() + labs(subtitle = "Metrópolis") +g_indep <- sims_indep_tbl %>% + filter(iter_num < 500) %>% + ggplot(aes(x = iter_num, y = theta)) + + geom_line() + labs(subtitle = "Independientes") +g_metropolis + g_indep +``` + +Donde vemos claramente que las simulaciones de metropolis están autocorrelacionadas: +la siguiente simulación depende de la anterior. Esto define una cadena de Markov. + +En cualquiera de los dos casos, como vimos en los histogramas de arriba, +las simulaciones "visitan" cada parte [0,1] de manera proporcional a la densidad, +de manera que podemos usar ambos tipos de simulaciones para aproximar la integral +o cantidad que nos interesa. Por ejemplo, la media posterior es: + +```{r} +media_1 <- sims_metro_tbl %>% summarise(media_post = mean(theta)) %>% pull(media_post) +media_2 <- sims_indep_tbl %>% summarise(media_post = mean(theta)) %>% pull(media_post) +media_exacta <- 3/(3 + 2) +tibble(metodo = c("sim Metrópolis", "sim Independiente", "exacto"), + media_post = c(media_1, media_2, media_exacta)) +``` + + + +```{block, type='mathblock'} +Supongamos que queremos simular de una distribución $p(\theta)$, +pero sólo conocemos $p(\theta)$ módulo una constante. Bajo ciertas condiciones de regularidad: + +El **algoritmo Metrópolis** para la distribución $p(\theta)$ define una **cadena +de Markov** cuya distribución a largo plazo es $p(\theta)$. Esto implica que si +$\theta^{(1)},\theta^{(2)}, \ldots, \theta^{(M)}$ es una simulación de esta +cadena, y $M$ es suficientemente grande + +1. La distribución de las $\theta^{(i)}$ es aproximadamente $p(\theta)$, +2. Tenemos que +$$ \frac1M \sum_{m = 1}^M h(\theta^{(m)}) \to \int h(\theta)p(\theta)\, d\theta$$ +cuando $M\to \infty$ +``` + +**Observaciones**: + +1. Aunque hay distintas *condiciones de regularidad* que pueden funcionar, +generalmente el supuesto es que la cadena de Markov construída es +[ergódica](https://en.wikipedia.org/wiki/Ergodic_theory), y hay varias +condiciones que garantizan esta propiedad. Una condición simple, por ejemplo, es +que el soporte de la distribución $p(\theta)$ es un conjunto conexo del espacio +de parámetros. + +2. Más crucialmente, este resultado no dice qué tan grande debe ser $M$ para que +la aproximación sea buena. Esto depende de cómo es $p(\theta)$, y de la +distribución que se utiliza para obtener los saltos propuestos. Dependiendo de +estos dos factores, la convergencia puede ser rápida (exponencial) o tan lenta +que es infactible usarla. Más adelante veremos diagnósticos para descartar los +peores casos de falta de convergencia. + +## Ajustando el tamaño de salto {-} + +En el algoritmo Metrópolis, generalmente es importante escoger la +dispersión de la distribución que genera propuestas con cuidado. + +- Si la dispersión de la propuesta es demasiado grande, tenderemos a rechazar mucho, +y la convergencia será lenta. +- Si la dispersión de la propuesta es demasiado chica, tardaremos mucho tiempo +en explorar las distintas partes de la distribución objetivo. + +### Ejemplo {-} + +Supongamos que queremos simular usando metróplis de una distribución +$\textsf{Gamma}(20, 100)$. Abajo vemos la forma de esta distribución: + +```{r} +sim_indep <- tibble(theta = rgamma(10000, 20, 100)) +ggplot(sim_indep, aes(x = theta)) + geom_histogram() +``` + +```{r} +# logaritmo de densidad no normalizada +log_f_dist <- function(x) 210 + dgamma(x, 20, 100, log = TRUE) +# iterar +iterador_metro_chico <- crear_metropolis(log_f_dist, sigma_salto = 0.001) +sims_chico_tbl <- iterador_metro_chico(c(theta = 0.02), 50000) +g_sim <- ggplot(sims_chico_tbl %>% filter(iter_num < 3000), aes(x = iter_num, y = theta)) + geom_line() + ylim(c(0, 0.5)) +dist_bplot <- ggplot(tibble(x = rgamma(10000, 20, 100)), aes(y = x, x = "a")) + geom_violin() + ylab("") + ylim(0, 0.5) +g_sim + dist_bplot + plot_layout(widths = c(5, 1)) +``` + +Nótese que después de 5 mil iteraciones estamos muy lejos de tener una muestra +que se aproxime a la distribución objetivo. Empezamos en un lugar bajo, y la +cadena sólo ha ido lentamente hacia las zonas de alta densidad. *Cualquier +resumen con esta cadena estaría fuertemente sesgado* al valor donde iniciamos la +iteración. Decimos que la cadena todavía no *mezcla* en las primeras 5 mil +iteraciones. + +Ahora vemos qué pasa si ponemos el tamaño de salto demasiado grande: + +```{r} +set.seed(831) +iterador_metro_grande <- crear_metropolis(log_f_dist, sigma_salto = 20) +sims_grande_tbl <- iterador_metro_grande(c(theta = 0.02), 50000) +g_sim <- ggplot(sims_grande_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)) +``` + +En este caso, la cadena se *atora* muchas veces, pues las propuestas tienen +probabilidad muy baja, y tendemos a tener una tasa de rechazos muy alta. Esto +quiere decir que la información que tenemos acerca de la posterior es +relativamente poca, pues muchos datos son repeticiones del mismo valor. +*Cualquier resumen con esta cadena podría estar muy lejos del verdadero valor,* +pues su varianza es alta - otra corrida se "atoraría" en otros valores +distintos. + +Nótese que cualquiera de estas cadenas, si la corremos suficientemente tiempo, +nos daría resultados buenos. Sin embargo, el número de simulaciones puede ser +infactible. + +Un valor intermedio nos dará mucho mejores resultados: + +```{r} +set.seed(831) +iterador_metro_apropiada <- crear_metropolis(log_f_dist, sigma_salto = 0.1) +sims_tbl <-iterador_metro_apropiada(c(theta = 0.02), 50000) +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. + +Comparemos cómo saldría por ejemplo la media posterior aproximada según +los tres métodos: + +```{r} +estimaciones_media <- map_dfr( + list(sims_chico_tbl, sims_grande_tbl, sims_tbl), + ~ filter(.x, iter_num < 3000) %>% + summarise(media = mean(theta))) %>% + mutate(tipo = c("salto chico", "salto grande", "salto apropiado")) +estimaciones_media %>% bind_rows(tibble(tipo = "exacta", media = 20/100)) %>% + select(tipo, media) +``` + +Veamos otra corrida: + +```{r} +set.seed(6222131) +sims_chica_tbl <- iterador_metro_chico(c(theta = 0.02), 5000) +sims_grande_tbl <- iterador_metro_grande(c(theta = 0.02), 5000) +estimaciones_media <- map_dfr( + list(sims_chica_tbl, sims_grande_tbl, sims_tbl), + ~ filter(.x, iter_num < 3000) %>% + summarise(media = mean(theta))) %>% + mutate(tipo = c("salto chico", "salto grande", "salto apropiado")) +estimaciones_media %>% bind_rows(tibble(tipo = "exacta", media = 20/100)) %>% + select(tipo, media) +``` + +```{block, type='ejercicio'} +Repite este proceso varias veces. Verifica que: + + - Si el tamaño de paso es muy chico, las estimaciones de la media tienen sesgo alto. + - Si el tamaño de paso es muy grande, las estimaciones tienen varianza alta. + - Si el tamaño de paso es adecuado, obtenemos buena precisión en la estimación de la media posterior. + - Explica estos tres casos en términos de la convergencia de las realizaciones + de la cadena de Markov. Explica cómo afecta a cada caso el valor inicial de las + simulaciones de Metrópolis. + - Repite para otra estadística, como la desviación estándar o el rangon intercuartil. +``` + + + +## Metrópolis con varios parámetros {-} + +Ahora aplicaremos el algoritmo Metrópolis cuando tenemos varios parámetros. La idea +es la misma, pero nuestra distribución de salto debe ser multivariada. Una selección +usual es usando saltos normales independientes para cada parámetro, es decir, la + normal multivariada con matriz de varianza y covarianza diagonal. + +### Ejemplo: el modelo normal {-} + +Veremos cómo simular con Metrópolis para el problema de los cantantes. +Sabemos como calcular la posterior: + +```{r} +crear_log_posterior_norm <- function(x = datos, m_0, n_0, a, b){ + # calcula log_posterior + log_posterior <- function(mu, sigma){ + log_verosim <- sum(dnorm(x, mu, sigma, log = TRUE)) + tau <- 1 / sigma^2 + log_inicial <- + dgamma(tau, a, b, log = TRUE) + + dnorm(mu, mu_0, sigma/sqrt(n_0), log = TRUE) + log_p <- log_verosim + log_inicial + log_p + } + log_posterior +} +``` + + +```{r} +# parametros de inicial y datos +a <- 3 +b <- 140 +mu_0 <- 175 +n_0 <- 5 +set.seed(3413) +cantantes <- lattice::singer %>% + mutate(estatura_cm = round(2.54 * height)) %>% + filter(str_detect(voice.part, "Tenor")) %>% + sample_n(20) +``` + +Vemos cómo se ven las primeras iteraciones de nuestra cadena de Markov: + +```{r} +log_p <- crear_log_posterior_norm(cantantes$estatura_cm, mu_0, n_0, a, b) +log_post <- function(pars) { log_p(pars[1], pars[2]) } +set.seed(823) +metro_normal <- crear_metropolis(log_post, sigma_salto = 0.5) +sim_tbl <- metro_normal(c(mu = 172, sigma = 3), 50000) +ggplot(sim_tbl %>% filter(iter_num < 100), + aes(x = mu, y = sigma)) + + geom_path() + + geom_point() +``` + +Y ahora vemos todas las simulaciones: + +```{r} +g_normal <- ggplot(sim_tbl, aes(x = mu, y = sigma)) + + geom_point(alpha = 0.05)+ coord_equal() + ylim(c(0, 14)) +g_normal +``` + +Y las medias posteriores son: + +```{r} +sim_tbl %>% summarise(across(is_double, mean)) +``` + +### Ejemplo: observaciones normales, no conjugado {-} + +Arriba repetimos el análisis conjugado usando Metrópolis. Aunque +ya no es necesario usar el modelo conjugado, y podemos poner +iniciales que sean más intuitivas y acorde con nuestro conocimiento +existente. + +Por ejemplo, podemos poner $p(\mu, \sigma) = p(\mu)p(\sigma)$, donde la densidad +de $\mu \sim \mathsf{N}(175, 2)$ y $\sigma \sim \mathsf{U}[2, 20].$ Igual que +antes, la verosimilitud $p(x|\mu, \sigma)$ es normal con media $\mu$ y +desviación estándar $\sigma.$ + +Escribimos la posterior: + +```{r} +crear_log_posterior <- function(x, m_0, sigma_0, inf, sup){ + # calcula log_posterior + log_posterior <- function(mu, sigma){ + log_verosim <- sum(dnorm(x, mu, sigma, log = TRUE)) + log_inicial <- + dunif(sigma, inf, sup, log = TRUE) + + dnorm(mu, mu_0, sigma_0, log = TRUE) + log_p <- log_verosim + log_inicial + log_p + } + log_posterior +} +``` + +```{r} +log_p <- crear_log_posterior(cantantes$estatura_cm, 175, 3, 2, 20) +log_post <- function(pars) { log_p(pars[1], pars[2]) } +``` + + +```{r, fig.width = 7} +set.seed(8231) +metro_normal <- crear_metropolis(log_post, sigma_salto = 0.5) +sim_tbl <- metro_normal(c(mu = 172, sigma = 5), 50000) +g_normal_2 <- ggplot(sim_tbl, aes(x = mu, y = sigma)) + + geom_point(alpha = 0.05) + coord_equal() + ylim(c(0, 14)) +g_normal + g_normal_2 +``` +Los resultados son similares, pero en +nuestras estimaciones bajo el segundo modelo, la $\sigma$ está +concentrada en valores un poco más bajos que el modelo normal-gamma inversa. +Las medias posteriores son: + +```{r} +sim_tbl %>% summarise(across(is.numeric, mean)) +``` + +Nótese que la inicial para el modelo normal-gamma inversa pone muy poca +probabilidad para valores bajos de $\sigma$, mientras que el segundo +modelo hay un 10% de probabilidad de que la $\sigma$ sea menor que 4. + +```{r} +tau <- rgamma(5000, 3, 150) +sigma <- 1/sqrt(tau) +quantile(sigma, c(0.01,0.1, 0.9, 0.99)) +quantile(runif(5000, 2, 25), c(0.01,0.1, 0.9, 0.99)) +``` + + + +### Ejemplo: exámenes {-} + +Recordamos un ejemplo que vimos en la sección de máxima verosimilitud. +Supongamos que en una población de estudiantes tenemos dos tipos: unos llenaron +un examen de opción múltiple al azar (1 de 5), y otros contestaron las +preguntas intentando sacar una buena calificación. Suponemos que una vez que +conocemos el tipo de estudiante, todas las preguntas tienen la misma +probabilidad de ser contestadas correctamente, de manera independiente. El +modelo teórico está representado por la siguiente simulación: + +```{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 +} +``` + +Y una muestra se ve como sigue: + +```{r} +set.seed(12) +muestra <- map_dbl(1:200, ~ sim_formas(0.35, 0.5)) +qplot(muestra) +``` + +Supongamos que no conocemos la probabildad de contestar correctamente ni la +proporción de estudiantes que contestó al azar. ¿Como estimamos estas dos cantidades? + +La verosimilitud la escribimos en el ejercicio anterior en la sección de máxima verosimilitud, está dada, para las repuestas de un estudiante, por: + +$$p(X = k|\theta_{azar}, \theta_{corr}) \propto \theta_{azar}(1/5)^k(4/5)^{10-k} + +(1-\theta_{azar})\theta_{corr}^k(1-\theta_{corr})^{10-k}$$ + +Suponiendo que todas las preguntas tienen la misma dificultad, y que +los estudiantes que estudiaron son homogéneos (podemos discutir qué haríamos +para introducir heterogeneidad que típicamente observaríamos). + +Creemos que la mayoría de los estudiantes no contesta al azar, así que pondremos +como inicial + +$$\theta_{azar} \sim \mathsf{Beta}(1, 5)$$ + +```{r} +qbeta(c(0.1, 0.9), 1, 5) %>% round(2) +``` +Ahora tenemos que pensar en la probabilidad $\theta_{corr}$ para los estudiantes +que sí estudiaron. Imaginemos que lo probamos con un estudiante que +sabemos que sí estudió, y obtuvo un porcentaje de correctos de 7/10, Podríamos +poner entonces (vimos 10 intentos, con 3 fracasos y 7 éxitos): + +$$\theta_{corr} \sim \mathsf{Beta}(7, 3)$$ +Finalmente, necesitamos la conjunta inicial. Pondremos +$$p(\theta_{azar},\theta_{corr}) = p(\theta_{azar})p(\theta_{corr})$$ +con lo que expresamos que inicialmente no creemos que estos dos parámetros estén +relacionados. Si pensáramos, por ejemplo, que cuando hacemos exámenes difíciles +menos estudiantes estudian, entonces deberíamos intentar otra conjunta. + +Escribimos el producto de la verosimilitud con la inicial: + +```{r} +crear_log_posterior <- function(x){ + + log_posterior <- function(theta_azar, theta_corr){ + log_verosim <- sum(log(theta_azar * dbinom(x, 10, 1/5) + + (1 - theta_azar) * dbinom(x, 10, theta_corr))) + log_inicial <- dbeta(theta_azar, 1, 5, log = TRUE) + + dbeta(theta_corr, 7, 3, log = TRUE) + log_post <- log_verosim + log_inicial + log_post + } + log_posterior + +} +``` + +Creamos la función de verosimilitud con los datos + +```{r} +log_p <- crear_log_posterior(muestra) +log_post <- function(pars) { log_p(pars[1], pars[2]) } +set.seed(8231) +metro_examenes <- crear_metropolis(log_post, sigma_salto = 0.02) +sim_tbl <- metro_examenes(c(theta_azar = 0.5, theta_corr = 0.5), 20000) +g_1 <- ggplot(sim_tbl, aes(x = theta_azar, y = theta_corr)) + + geom_point(alpha = 0.05) + coord_equal() +g_1 +``` +Nótese que hay cierta correlación entre las dos proporciones, y esto produce +intervalos posteriores relativamente amplios. Esto es de esperarse, pues +los datos son consistentes con una proporción relativamente chica de +estudiantes que contestan al azar, y tasas de correctos más altas entre los +que sí estudian, y una proporción más grande de respuestas al azar con +tasas de correctos más altas. + +```{r} +f <- c(0.05, 0.5, 0.95) +sim_tbl %>% + pivot_longer(-iter_num, names_to = "parametro", values_to = "valor") %>% + group_by(parametro) %>% + summarise(cuantil = quantile(valor, f), f = f) %>% + mutate(cuantil = round(cuantil, 2)) %>% + pivot_wider(names_from = f, values_from = cuantil) +``` + + +## Muestreador de Gibbs {-} + +El algoritmo de Metrópolis es muy general y se puede aplicar a una gran variedad +de problemas. Sin embargo, afinar los parámetros de la distribución propuesta +para que el algoritmo funcione correctamente puede ser complicado. El +muestredor de Gibbs no necesita de una distribución propuesta y por lo tanto +no requiere afinar estos parámetros. + +**Para implementar un muestreador de Gibbs se necesita ser capaz de generar +muestras de la distribución posterior condicional a cada uno de los +parámetros individuales.** Esto es, el muestreador de Gibbs permite generar +muestras de la posterior: +$$p(\theta_1,...,\theta_p|x)$$ +siempre y cuando podamos generar valores de todas las distribuciones +condicionales: +$$\theta_k \sim p(\theta_k|\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p,x).$$ + +El proceso del muestreador de Gibbs es una caminata aleatoria a lo largo del +espacio de parámetros. La caminata inicia en un punto arbitrario y en cada +tiempo el siguiente paso depende únicamente de la posición actual. Por tanto +el muestredor de Gibbs es un proceso cadena de Markov vía Monte Carlo. La +diferencia entre Gibbs y Metrópolis radica en como se deciden los pasos. + + +**Muestreador Gibbs** + +En cada punto de la caminata se selecciona uno de los +componentes del vector de parámetros (típicamente se cicla en orden): + +1. Supongamos que se selecciona el parámetro $k$-ésimo después de haber +modificado los $k-1$ anteriores, entonces obtenemos un nuevo valor para este +parámetro generando una simulación de la distribución condicional +$$\theta_k^{(i+1)} \sim p(\theta_k|\theta_1^{(i+1)},\ldots,\theta_{k-1}^{(i+1)},\theta_{k+1}^{(i)},\ldots,\theta_p^{(i)},x)$$ + +2. El nuevo valor $\theta_k^{(i+1)}$ junto con los valores +$\theta_1^{(i+1)},\ldots,\theta_{k-1}^{(i+1)},\theta_{k+1}^{(i)},\ldots,\theta_p^{(i)}$ +constituyen la nueva posición en la caminata aleatoria. + +3. Seleccionamos una nueva componente $\theta_{k+1}^{(i+1)}$ y repetimos el proceso. + + +El muestreador de Gibbs es útil cuando no podemos determinar de manera analítica +la distribución conjunta y no se puede simular directamente de ella, pero sí +podemos determinar todas las distribuciones condicionales y simular de ellas. + +### Ejemplo: dos proporciones {-} + +Supongamos que queremos evaluar el balanceo de *dos dados* de 20 lados que +produce una fábrica. En particular, evaluar la probabilidad de tirar un 20, y +quizá escoger el dado que nos de mayor probabilidad de tirar un 20. + +Tiramos cada dado $n$ veces, y denotamos por $X_1$ y $X_2$ el número +de 20's que tiramos en cada ocasión. El modelo de datos está dado por +$$p(x_1, x_2|\theta_1, \theta_2)\propto \theta_1^{x_1}(1-\theta_1)^{n - x_1}\theta_2^{x_2}(1-\theta_2)^{n - x_2},$$ +que es el producto de dos densidades binomiales, pues suponemos que +las tiradas son independientes cuando conocemos los parámetros +$\theta_1$ y $\theta_2$. + +Ahora ponemos una inicial +$$p(\theta_i)\sim \mathsf{Beta}(100, 1900)$$ + +y aquí están las razones de nuestra elección: +```{r} +media <- 1/20 +k <- 2000 +a <- media * k +b <- (1 - media) * k +c(a,b) +qbeta(c(0.05, 0.95), a, b) %>% round(3) +``` +y suponemos que + +$$p(\theta_1,\theta_2) = p (\theta_1)p(\theta_2)$$ + +es decir, apriori saber el desempeño de un dado no nos da información adicional del otro (esto podría no ser cierto, por ejemplo, si el defecto es +provocado por la impresión del número 20). + +Por lo tanto, la posterior es + +$$p(\theta_1,\theta_2|x_1, x_2)\propto \theta_1^{x_1+100-1}(1-\theta_1)^{n - x_1 + 1900-1}\theta_2^{x_2+100 -1}(1-\theta_2)^{n - x_2 + 1900-1}$$ + +Ahora consideramoso qué pasa cuando conocemos $\theta_2$ y los datos. Pensamos en todo lo que no sea $\theta_1$ como constante de modo que nos +queda: + +$$p(\theta_1 | \theta_2, x) \propto \theta_1^{x_1+100 -1}(1-\theta_1)^{n - x_1 + 1900 -1}$$ +que es $\mathsf{Beta}(x_1 + 100, n - x_1 + 1900)$, y por la misma razón, + +$$p(\theta_2 | \theta_1, x) \propto \theta_2^{x_2+100-1}(1-\theta_2)^{n - x_2 + 1900-1}$$ + +que también es es $\mathsf{Beta}(x_1 + 100, n - x_1 + 1900)$ + +De hecho, estas condicionales son fáciles de deducir de otra manera: en +realidad estamos haciendo dos experimentos separados (pues suponemos +que las iniciales son independientes y las pruebas también), así que +podriamos usar el análisis Beta-Binomial para cada uno de ellos. En realidad no es necesario usar MCMC para este ejemplo. + +Usaremos esta función para hacer nuestras iteraciones de Gibbs: + + +```{r} +iterar_gibbs <- function(pasos, n, x_1, x_2){ + iteraciones <- matrix(0, nrow = pasos + 1, ncol = 2) # vector guardará las simulaciones + iteraciones[1, 1] <- 0.5 # valor inicial media + colnames(iteraciones) <- c("theta_1", "theta_2") + # Generamos la caminata aleatoria + for (j in seq(2, pasos, 2)) { + # theta_1 + a <- x_2 + 100 - 1 + b <- n - x_2 + 1900 - 1 + iteraciones[j, "theta_2"] <- rbeta(1, a, b) # Actualizar theta_1 + iteraciones[j, "theta_1"] <- iteraciones[j-1, "theta_1"] + # theta_2 + a <- x_1 + 100 - 1 + b <- n - x_1 + 1900 - 1 + iteraciones[j + 1, "theta_1"] <- rbeta(1, a, b) # Actualizar theta_1 + iteraciones[j + 1, "theta_2"] <- iteraciones[j, "theta_2"] + } + iteraciones +} +``` + + +Y supongamos que estamos comparando los dados de dos compañías: Chessex y GameScience. Tiramos cada dado 10 mil veces, y obtenemos: + +```{r} +# Datos de https://www.awesomedice.com/blogs/news/d20-dice-randomness-test-chessex-vs-gamescience +n <- 10000 +x_1 <- 408 # Chessex, alrededor de 0.85 dólares por dado +x_2 <- 474 # GameScience, alrededor 1.60 dólares por dado +``` + +E iteramos: + +```{r} +iteraciones <- iterar_gibbs(20000, n, x_1, x_2) %>% + as_tibble() %>% + mutate(iter_num = row_number()) +head(iteraciones) +``` + +```{r} +ggplot(filter(iteraciones, iter_num > 1000, iter_num< 1050), + aes(x = theta_1, y = theta_2)) + + geom_path(alpha = 0.3) + geom_point() +``` + +```{r} +g_1 <- ggplot(iteraciones, + aes(x = theta_1, y = theta_2)) + + geom_path(alpha = 0.3) + geom_point() +g_2 <- ggplot(iteraciones %>% filter(iter_num > 10), + aes(x = theta_1, y = theta_2)) + + geom_path(alpha = 0.3) + geom_point() + + geom_abline(colour = "red") + + geom_point(data= tibble(theta_1=1/20, theta_2=1/20), colour = "red", size = 5) +g_1 + g_2 +``` + +Notamos el dado de Cheesex no es consistente con 1/20 de tiros de 20s, +pero el dado de GameScience sí lo es. De este gráfica vemos +que Cheesex está sesgado hacia abajo, así que deberíamos escoger +el dado de GameScience + +Podemos ver directamente cómo se distribuye la diferencia $\theta_1 - \theta_2$. Cualquier estadística es fácil de evaluar, pues simplemente +la calculamos para cada simulación y después resumimos: + +```{r} +iteraciones <- iteraciones %>% + mutate(dif = theta_1 - theta_2) +ggplot(iteraciones %>% filter(iter_num > 10), aes(x = dif)) + + geom_histogram(bins = 100) + + geom_vline(xintercept = 0, colour = "red") +``` + +Y vemos que es altamente probable que el dado de Cheesex produce +más 20's que el dado de GameScience. + +```{r} +iteraciones %>% mutate(theta_1_mayor = dif > 0) %>% + summarise(prob_theta_1_mayor = mean(theta_1_mayor)) +``` + + +Finalmente, verificamos nuestro modelo y cuánto aprendimos. Podemos +hacerlo simulando de la inicial y comparando con la posterior: + +```{r} +inicial_tbl <- tibble(theta_1 = rbeta(20000, 100, 1900), + theta_2 = rbeta(20000, 100, 1900), + dist = "inicial") +posterior_tbl <- iteraciones %>% filter(iter_num > 10) %>% + mutate(dist = "posterior") +sims_tbl <- bind_rows(inicial_tbl, posterior_tbl) +ggplot(sims_tbl, aes(x = theta_1, y = theta_2, colour = dist)) + + geom_point(alpha = 0.2) +``` +donde vemos que el resultado que obtuvimos es razonablemente consistente +con nuestra información inicial, y las 10 mil tiradas de dado +fueron altamente informativas. + +```{block2, type="ejercicio"} + +- ¿Qué crees que pasaría si sólo hubieramos tirado 40 veces cada +dado? ¿Qué tanto habríamos aprendido? Puedes usar datos +simulados y repetir este ejercicio. + +- Puedes examinar los resultados para cada cara con los datos originales. +Un modelo apropiado es el Dirichlet-Multinomial. + +``` + + + + + +### Ejemplo: Modelo normal no conjugado {-} + +Retomemos el caso de observaciones normales, supongamos que tenemos una muestra +$X_1,...,X_n$ de observaciones independientes e identicamente distribuidas, +con $X_i \sim \mathsf{N}(\mu, \sigma^2)$. + +Usaremos iniciales distintas al modelo anterior: + +$$p(\mu, \sigma^2) = p(\sigma^2)p(\mu)$$ +con $\mu$ $\mathsf{N}(\mu_0, \sigma_0)$ y $\tau = 1/\sigma^2$ con distribución +$\mathsf{Gamma}(a,b)$. Esto no nos da el modelo conjugado que vimos antes (nota +la diferencia de la especificación de la inicial conjunta). + +Comenzamos por escribir + +$$p(\mu, \sigma^2|x) \propto \frac{1}{{\sigma^{n/2}}} \exp(-\sum\frac{(x_i-\mu)²}{2\sigma^2}) \exp(- \frac{(\mu - \mu_0)^2}{2\sigma_0^2}) \frac{1}{(\sigma^2)^{a + 1}}\exp (-\beta/\sigma^2 )$$ + +Comenzamos analizando $p(\mu|\sigma^2, x)$. Por la ecuación de arriba, e ignorando +los términos que **no** dependen de $\mu$: +$$p(\mu|\sigma^2, x) \propto \exp [ - \sum_i (\frac{(\mu - x_i)^2}{2\sigma^2} - \frac{(\mu - \mu_0)^2}{2n\sigma_0^2})]$$ +que es una distribución normal (completa cuadrados): + +$$\mu|\sigma^2,x \sim \mathsf{N}\bigg(\frac{\sigma^2}{\sigma^2 + n\sigma_0^2}\mu_0 + \frac{n\sigma_0^2}{\sigma^2 + n \sigma_0^2}\bar{x}, \frac{\sigma \sigma_0}{\sqrt{\sigma^2 + n\sigma_0^2}}\bigg)$$ + +Ahora consideramos $p(\sigma^2|mu,x)$. Ignoramos en $p(\mu,\sigma^2|x)$ los términos que *no** dependen de $\sigma^2$: + +$$p(\sigma^2|\mu, x) \propto \frac{1}{\sigma^{n/2}} \exp(-\sum\frac{(x_i-\mu)²}{2\sigma^2}) \frac{1}{(\sigma^2)^{a + 1}}\exp (-\beta/\sigma^2)$$ +que simplificando da + +$$ = \frac{1}{\sigma^{n/2 + a + 1}}\exp( -\frac{\beta +\frac{1}{2}\sum(x_i - \mu)^2}{\sigma^2} )$$ +de modo que + +$$\sigma^2|\mu, x \sim \mathsf{GI}(a +n/2, b + \frac{1}{2}\sum(x_i -\mu)^2)$$ + +### Ejemplo {-} + +Usaremos este muestreador para el problema de la estaturas de +los tenores. Comenzamos definiendo las distribuciones iniciales: + +* $\mu \sim \mathsf{N}(175, 3)$ + +* $\tau = 1/\sigma^2 \sim \mathsf{GI}(3, 150)$, esto es $a = 3$ y $b = 150$. + +Escribimos el muestreador de Gibbs. + +```{r} +n <- 20 +x <- cantantes$estatura_cm +m <- 175; sigma_0 <- 3; alpha <- 3; beta <- 150 # parámetros de iniciales +pasos <- 20000 +iteraciones <- matrix(0, nrow = pasos + 1, ncol = 2) # vector guardará las simulaciones +iteraciones[1, 1] <- 0 # valor inicial media +colnames(iteraciones) <- c("mu", "sigma") +# Generamos la caminata aleatoria +for (j in seq(2, pasos, 2)) { + # sigma^2 + mu <- iteraciones[j - 1, "mu"] + a <- n / 2 + alpha + b <- sum((x - mu) ^ 2) / 2 + beta + iteraciones[j, "sigma"] <- sqrt(1/rgamma(1, a, b)) # Actualizar sigma + iteraciones[j, "mu"] <- iteraciones[j-1, "mu"] + # mu + sigma <- iteraciones[j, "sigma"] + media <- (n * sigma_0^2 * mean(x) + sigma^2 * m) / (n * sigma_0^2 + sigma^2) + varianza <- sigma^2 * sigma_0^2 / (n * sigma_0^2 + sigma^2) + iteraciones[j+1, "mu"] <- rnorm(1, media, sd = sqrt(varianza)) # actualizar mu + iteraciones[j+1, "sigma"] <- iteraciones[j, "sigma"] +} +caminata <- data.frame(pasos = 1:pasos, mu = iteraciones[1:pasos, "mu"], + sigma = iteraciones[1:pasos, "sigma"]) +caminata_g <- caminata %>% + gather(parametro, val, mu, sigma) %>% + arrange(pasos) +``` + + +Veamos primero algunos pasos: + +```{r} +ggplot(filter(caminata, pasos > 1000, pasos< 1010), + aes(x = mu, y = sigma)) + + geom_path(alpha = 0.3) + geom_point() +``` + +Donde vemos cómo en cada iteración se actualiza un solo parámetro. +Una alternativa es conservar únicamente ciclos completos de la +caminata u esto es lo que hacen varios programas que implementan Gibbs, sin +embargo ambas cadenas (cadenas completas y conservando únicamente ciclos +completos) convergen a la misma distribución posterior. + +Si tomamos iteraciones completas: + +```{r} +ggplot(filter(caminata, pasos > 1000, pasos< 1020, pasos %% 2 == 0), + aes(x = mu, y = sigma)) + + geom_path(alpha = 0.3) + geom_point() +``` +Y ahora vemos cómo se ven las simulaciones: + +```{r} +ggplot(filter(caminata, pasos > 1000, pasos< 10000, pasos %% 2 == 0), + aes(x = mu, y = sigma)) + geom_point(alpha = 0.1) +``` + +Y el diagnóstico de cada cadena: + +```{r} +ggplot(filter(caminata_g, pasos > 15000), aes(x = pasos, y = val)) + + geom_path(alpha = 0.3) + + facet_wrap(~parametro, ncol = 1, scales = "free") + + scale_y_continuous("") +``` +Estas cadenas parecen estar mezclando bien. Podemos resumirlas: + +```{r, fig.width = 4} +ggplot(filter(caminata_g, pasos > 5000), aes(x = val)) + + geom_histogram(fill = "gray") + + facet_wrap(~parametro, ncol = 1, scales = "free") +``` +```{r} +caminata_g %>% + filter(pasos > 1000) %>% # eliminamos la etapa de calentamiento + group_by(parametro) %>% + summarise( + mean(val), + sd(val), + median(val) + ) %>% + mutate(across(is_double, round, 2)) +``` +Y obtenemos un resultado similar a los anteriores. + + + +## Conclusiones y observaciones Metrópolis y Gibbs {-} + +* Una generalización del algoritmo de Metrópolis es Metrópolis-Hastings. + + El algoritmo de Metrópolis es como sigue: + 1. Generamos un punto inicial tal que $p(\theta)>0$. + 2. Para $i = 1,2,...$ + + Se propone un nuevo valor $\theta^*$ con una distribución + propuesta $g(\theta^*|\theta^{(i)})$ es común que $g(\theta^*|\theta^{(i)})$ sea una normal centrada en + $\theta^{(i)}$. + 3. Calculamos la probabilidad de aceptación + + $$\alpha=\min\bigg\{\frac{p(\theta^*)}{p(\theta^{(i)})},1\bigg\},$$ + y aceptamos $\theta^*$ con probabilidad $p_{mover}$. Es así que el + algorito requiere que podamos calcular el cociente en $p_{mover}$ para todo + $\theta^{(i)}$ y $\theta^*$, así como simular de la distribución propuesta + $g(\theta^*|\theta^{(i)})$, adicionalmente debemos poder generar valores + uniformes para decidir si aceptar/rechazar. + + En el caso de **Metrópolis** un requerimiento adicional es que la distribución + propuesta $g(\theta_{a}|\theta_b)$ debe ser simétrica, es decir + $g(\theta_{a}|\theta_b) = g(\theta_{b}|\theta_a)$ para todo $\theta_{a}$, + $\theta_{b}$. + + **Metrópolis-Hastings** generaliza Metrópolis, eliminando la restricción de + simetría en la distribución propuesta $g(\theta_{a}|\theta_b)$, sin embargo para + corregir por esta asimetría debemos calcular $\alpha$ como sigue: + + $$\alpha=\min\bigg\{ \frac{p(\theta^*)}{g(\theta^*|\theta^{(i)})} \cdot \frac{g(\theta^{(i)}|\theta^*)}{p(\theta^{(i)})},1\bigg\}$$ +La generalización de Metrópolis-Hastings puede resultar en algoritmos más +veloces. + +* Se puede ver Gibbs como una generalización de Metrópolis-Hastings, cuando +estamos actualizando un componente de los parámetros, la distribución propuesta +es la distribución posterior para ese parámetro, por tanto siempre es aceptado. + +* Comparado con Metrópolis, Gibbs tiene la ventaja de que no se necesita afinar +los parámetros de una distribución propuesta (o seleccionar siquiera una +distribución propuesta). Además que no hay pérdida de simulaciones debido a +rechazo. Por su parte, la desventaja debemos conocer las distribuciones +condicionales y poder simular de ellas. + +* En el caso de modelos complicados se utilizan combinaciones de Gibbs y +Metrópolis. Cuando se consideran estos dos algoritmos Gibbs es un método más +simple y es la primera opción para modelos condicionalmente conjugados. Sí solo +podemos simular de un subconjunto de las distribuciones condicionales +posteriores, entonces podemos usar Gibbs siempre que se pueda y Metrópolis +unidimensional para el resto, o de manera más general separamos en bloques, un +bloque se actualiza con Gibbs y otro con Metrópolis. + +* El algoritmo de Gibbs puede *atorarse* cuando hay correlación alta entre los +parámetros, reparametrizar puede ayudar, o se pueden usar otros algoritmos. + +* [JAGS](http://mcmc-jags.sourceforge.net) (Just Another Gibbs Sampler), WinBUGS +y OpenBUGS son programas que implementan métodos MCMC para generar simulaciones +de distribuciones posteriores. Los paquetes `rjags` y `R2jags` permiten ajustar +modelos en JAGS desde `R`. Es muy fácil utilizar estos programas pues uno +simplemente debe especificar las distribuciones iniciales, la verosimilitud y +los datos observados. Para aprender a usar JAGS se puede revisar la sección +correspondiente en las [notas de 2018](https://tereom.github.io/est-computacional-2018/jags.html), +ahora nos concentraremos en el uso de Stan. + +## HMC y Stan {-} + +> It appears to be quite a general principle that, whenever there is a +randomized way of doing something, then there is a nonrandomized way that +delivers better performance but requires more thought. -E.T. Jaynes + +`Stan` es un programa para generar muestras de una distribución posterior de los +parámetros de un modelo, el nombre del programa hace referencia a [Stanislaw Ulam (1904-1984)](https://en.wikipedia.org/wiki/Stanislaw_Ulam) que fue pionero en +los métodos de Monte Carlo. A diferencia de JAGS y BUGS, los pasos de la cadena +de Markov se generan con un método llamado *Monte Carlo Hamiltoniano* (HMC). HMC +es computacionalmente más costoso que Metrópolis o Gibbs, sin embargo, sus +propuestas suelen ser más eficientes, y por consiguiente no necesita muestras +tan grandes. En particular cuando se ajustan modelos grandes y complejos (por +ejemplo, con variables con correlación alta) HMC supera a otros. + + + + + +## Diagnósticos generales para MCMC {-} + +Cuando generamos una muestra de la distribución posterior usando MCMC, sin +importar el método (Metrópolis, Gibbs, HMC), buscamos que: + +1. Los valores simulados sean representativos de la distribución posterior. Esto implica que no deben estar influenciados por el valor inicial (arbitrario) y deben explorar todo el rango de la posterior, con suficientes retornos para evaluar cuánta masa hay en cada región. + +2. Debemos tener suficientes simulaciones de tal manera que las +estimaciones sean precisas y estables. + +3. Queremos tener un método eficiente para generar las simulaciones. + + +En la práctica intentamos cumplir lo más posible estos objetivos, pues aunque en principio los métodos MCMC garantizan que una cadena infinitamente larga logrará una representación perfecta, siempre debemos tener un criterio para cortar la cadena y evaluar la calidad de las simulaciones. + +### Representatividad {-} + +**Burn-in e iteraciones iniciales**- En primer lugar, en muchas ocasiones las condiciones iniciales de las +cadenas están en partes del espacio de parámetros que son "atípicos" en +términos de la posterior. Así que es común quitar algunas observaciones +iniciales (iteraciones de *burn-in*) para minimizar su efecto en +resúmenes posteriores. + + +Por ejemplo, para el ejemplo de los cantantes, podemos +ver que las iteraciones iniciales tienen como función principal +llegar a las regiones de probabilidad posterior alta: + +```{r} +log_p <- crear_log_posterior_norm(cantantes$estatura_cm, mu_0, n_0, a, b) +log_post <- function(pars) { log_p(pars[1], pars[2]) } +set.seed(823) +metro_normal <- crear_metropolis(log_post, sigma_salto = 0.5) +sim_tbl <- metro_normal(c(mu = 162, sigma = 1), 5000) +ggplot(sim_tbl %>% filter(iter_num < 500), aes(x = mu, y = sigma)) + geom_path(alpha = 0.5) + geom_point(aes(colour = iter_num)) +``` + +De modo que puede ser buena idea eliminar las +primeras iteraciones. En teoría, no es necesario hacer esto si +hacemos suficientes iteraciones, pues la cadena va a terminar +en su estado estable explorando la posterior. En la práctica, y +con pocas iteraciones, puede ayudar un poco a mejorar la precisión +numérica de las cantidades que queramos calcular. + +```{r, fig.width=7} +sim_g <- sim_tbl %>% pivot_longer(-iter_num, + names_to = "parametro", + values_to = "valor") +todas <- ggplot(sim_g, aes(x = iter_num, y = valor)) + + geom_line(alpha = 0.5) + + facet_wrap(~ parametro, ncol = 1, scales = "free_y") + + labs(subtitle = "Todas las simulaciones") +sin_burnin <- + sim_g %>% filter(iter_num > 200) %>% + ggplot(aes(x = iter_num, y = valor)) + + geom_line(alpha = 0.5) + + facet_wrap(~ parametro, ncol = 1, scales = "free_y") + + labs(subtitle = "Quitando 200 de burn-in") +todas + sin_burnin +``` + + + + +**Convergencia a estado límite**. Para determinar la convergencia es conveniente +realizar **más de una cadena**: buscamos ver si realmente se ha olvidado el +estado inicial, si las distribuciones de cada cadena son consistentes unas con +otras, y revisar que algunas cadenas no hayan quedado *atoradas* en regiones +inusuales del espacio de parámetros. + +Inicializamos las cadenas con valores al azar en rangos +razonables (por ejemplo simulando de la inicial): + +```{r} +set.seed(8513) +valores_iniciales <- tibble(mu_0 = rnorm(4, 160, 20), + sigma_0 = runif(4, 0, 20), + cadena = 1:4) +sims_tbl <- valores_iniciales %>% + mutate(sims = map2(mu_0, sigma_0, + ~ metro_normal(c(mu = .x, sigma = .y), 300) )) %>% + unnest(sims) + +ggplot(sims_tbl, aes(x = iter_num, y = sigma, colour = factor(cadena))) + + geom_line() +``` + +Y este es un ejemplo donde claramente las cadenas **no** han alcanzado +un estado estable: tienen muy distintas medias y varianzas. Por ejemplo: + + +```{r} +set.seed(83243) +sims_tbl <- valores_iniciales %>% + mutate(sims = map2(mu_0, sigma_0, + ~ metro_normal(c(mu = .x, sigma = .y), 20000) )) %>% + unnest(sims) + +ggplot(sims_tbl, aes(x = iter_num, y = sigma, colour = factor(cadena))) + + geom_line() +``` + +Y este resultado se ve mejor. La parte *transición* hacia las zonas +de alta probabilidad pasa antes de unas 1000 iteraciones. Podemos +hacer más simulaciones, o eliminar como *burn-in* las primiras iteraciones: + +```{r, fig.width = 6} +media_g <- ggplot(sims_tbl %>% filter(iter_num > 2000), + aes(x = iter_num, y = mu, colour = factor(cadena))) + + geom_line() +sigma_g <- ggplot(sims_tbl %>% filter(iter_num > 2000), + aes(x = iter_num, y = sigma, colour = factor(cadena))) + + geom_line() +media_g / sigma_g +``` + +Las gráficas anteriores nos ayudan a determinar si elegimos un periodo de +calentamiento adecuado o si alguna cadena está alejada del resto. + +Una vez que las cadenas están en estado estable, podemos usar +**todas** las simulaciones juntas para resumir: + +```{r} +head(sims_tbl) +# medias posteriores +sims_tbl %>% + summarise(mu = mean(mu), sigma = mean(sigma)) +``` + + + + + +Además de realizar gráficas podemos usar la medida de convergencia $\hat{R}$. La medida $\hat{R}$ se conoce como el **factor de reducción potencial de +escala** o *diagnóstico de convergencia de Gelman-Rubin*, esta es una estimación +de la posible reducción en la longitud de un intervalo de confianza si las +simulaciones continuaran infinitamente. $\hat{R}$ es aproximadamente la raíz +cuadrada de la varianza de todas las +cadenas juntas dividida entre la varianza dentro de cada cadena. Si $\hat{R}$ es +mucho mayor a 1 esto indica que las cadenas no se han mezclado bien. Una regla +usual es iterar hasta alcanzar un valor $\hat{R} \leq 1.1$ para todos los +parámetros. + +$$\hat{R} \approx \sqrt{\frac{\hat{V}}{W}}$$ + +donde $B$ es la varianza entre las cadenas, $W$ es la varianza dentro de las cadenas + +$$B = \frac{N}{M-1}\sum_m (\hat{\theta}_m - \hat{\theta})^2$$ +$$W = \frac{1}{M}\sum_m \hat{\sigma}_m^2$$ + +Y $\hat{V}$ es una estimación del varianza de posterior de $\theta$: + +$$\hat{V} = \frac{N-1}{N}W + \frac{M+1}{MN}B$$ +#### Ejemplo {-} + +En nuestro ejemplo anterior, tenemos + +```{r} +sims_tbl %>% + pivot_longer(mu:sigma, names_to = "parametro", values_to = "valor") %>% + group_by(parametro, cadena) %>% + summarise(media = mean(valor), num = n(), sigma2 = var(valor)) %>% + summarise(N = first(num), + M = n_distinct(cadena), + B = N * var(media), + W = mean(sigma2), + V_hat = ((N - 1) / N) * W + (M + 1)/(M * N) * B, + R_hat = sqrt(V_hat / W)) + +``` +Y verificamos que los valores de $\hat{R}$ son cercanos a uno, lo +cual indica que este diagnóstico es aceptable. Si hubiéramos +trabajado con las primeras 300 iteraciones + +```{r} +sims_tbl %>% + filter(iter_num < 300) %>% + pivot_longer(mu:sigma, names_to = "parametro", values_to = "valor") %>% + group_by(parametro, cadena) %>% + summarise(media = mean(valor), num = n(), sigma2 = var(valor)) %>% + summarise(N = first(num), + M = n_distinct(cadena), + B = N * var(media), + W = mean(sigma2), + V_hat = ((N - 1) / N) * W + (M + 1)/(M * N) * B, + R_hat = sqrt(V_hat / W)) + +``` +Y estos valores indican problemas en la convergencia de las cadenas. Es +necesario diagnosticar el problema, que en este caso resolvemos +incrementando el número de iteraciones. + + +### Precisión {-} + +Una vez que tenemos una muestra representativa de la +distribución posterior, nuestro objetivo es asegurarnos de que la muestra es lo suficientemente grande +para producir estimaciones estables y precisas de la distribución. + +Para ello usaremos el +**tamaño efectivo de muestra**, Si las simulaciones fueran independientes +$N_{eff}$ sería el número total de simulaciones; sin embargo, las simulaciones de MCMC suelen estar correlacionadas, de modo que cada iteración +de MCMC es menos informativa que si fueran independientes. + +**Ejemplo**: Si graficaramos simulaciones independientes, esperaríamos valores de +autocorrelación chicos: + +```{r} +acf(rgamma(1000,1,1)) +``` +Sin embargo, los valores que simulamos tienen el siguiente perfil de +autocorrelación: + +```{r} +sigma_metro_sims <- sims_tbl %>% filter(cadena==4) %>% pull(mu) +acf(sigma_metro_sims) +``` + +El tamaño efectivo de muestra nos dice qué tamaño de +muestra de observaciones independientes nos daría la misma información que las +simulaciones de la cadena. Una manera de manera relativamente simple de +estimarlo es: + +$$N_{eff} = \frac{N}{1+2\sum_{k=1}^\infty ACF(k)} $$ + +Usualmente nos gustaría obtener un tamaño efectivo de al menos $100$ (para +cálculo de medias y varianzas posteriores). Esta +cantidad usualmente se reporta en el software (con mejores estimaciones que +la de la fórmula de arriba), y es necesario checarlo. + +En nuestro ejemplo hacemos una aproximación como sigue: + +```{r} +calc_acf <- function(x){ + valores_acf <- acf(x, lag.max = 1000, plot = FALSE)$acf %>% as.numeric() + valores_acf[-1] +} +acf_tbl <- sims_tbl %>% + pivot_longer(mu:sigma, names_to = "parametro", values_to = "valor") %>% + group_by(parametro, cadena) %>% + summarise(N = n_distinct(iter_num), k = 1:1000, acf = calc_acf(valor)) %>% + summarise(N = first(N), N_eff = N / (1 + 2 * sum(acf))) +acf_tbl +``` + +Nótese que algunas cadenas tienen un tamaño efectivo de muestra relativamente +bajo para el número de iteraciones que hicimos. De cualquier forma, el agregado +sobre todas las cadenas es suficientemente grande para calcular resúmenes básicos: + +```{r} +acf_tbl %>% group_by(parametro) %>% + summarise(N = sum(N), N_eff = sum(N_eff)) +``` +Sin embargo, podemos hacer más simulaciones si es necesario, por ejemplo +para aproximar de manera apropiada percentiles en las colas. + + + + +### Eficiencia {-} + +Hay varias maneras para mejorar la eficiencia de un proceso MCMC: + + +* Paralelizar, no disminuimos el número de pasos en las simulaciones pero +podemos disminuir el tiempo que tarda en correr. + +* Cambiar la parametrización del modelo o transformar los datos. + +* Adelgazar la muestra cuando tenemos problemas de uso de memoria, + +consiste en guardar únicamente los $k$-ésimos pasos de la cadena y resulta +en cadenas con menos autocorrelación . + +### Recomendaciones generales {-} + +@gelman-hill recomienda los siguientes pasos cuando uno esta simulando de la +posterior: + +1. Cuando definimos un modelo por primera vez establecemos un valor bajo para +el número de iteraciones. La razón es que la mayor parte de las veces los +modelos no funcionan a la primera por lo que sería pérdida de tiempo dejarlo +correr mucho tiempo antes de descubrir el problema. + +2. Si las simulaciones no han alcanzado convergencia aumentamos las iteraciones +a $500$ ó $1000$ de tal forma que las corridas tarden segundos o unos cuantos +minutos. + +3. Si tarda más que unos cuantos minutos (para problemas del tamaño que +veremos en la clase) y aún así no alcanza convergencia +entonces _juega_ un poco con el modelo (por ejemplo intenta transformaciones lineales), para JAGS Gelman +sugiere más técnicas para acelerar la convergencia en el +capitulo $19$ del libro +*Data Analysis Using Regression and Multilevel/Hierarchical models*. En el +caso de Stan veremos ejemplos de reparametrización, y se puede leer más en +la [guía](https://mc-stan.org/docs/2_21/stan-users-guide/reparameterization-section.html). + +4. Otra técnica conveniente cuando se trabaja con bases de datos grandes +(sobre todo en la parte exploratoria) es trabajar con un +subconjunto de los datos, quizá la mitad o una quinta parte. + diff --git a/renv.lock b/renv.lock index e384f62..813f1b9 100644 --- a/renv.lock +++ b/renv.lock @@ -735,6 +735,17 @@ ], "Hash": "a1066cbc05caee9a4bf6d90f194ff4da" }, + "expm": { + "Package": "expm", + "Version": "1.0-0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "methods" + ], + "Hash": "a44b2810f36c1cda5d52eee6ec96cafa" + }, "fansi": { "Package": "fansi", "Version": "1.0.6",