From fe993338fb45958bfa8aa3d46cce4421bce2f958 Mon Sep 17 00:00:00 2001 From: Felipe Gonzalez Date: Tue, 12 Mar 2024 19:58:22 -0600 Subject: [PATCH] =?UTF-8?q?Agregar=20principio=20de=20diagn=C3=B3sticos?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- datos/wines_2012.csv | 181 +++++++++++++++++++++++++++++++++++++++++++ notas/08-mcmc.qmd | 176 ++++++++++++++++++++++++++++++++++++++++- renv.lock | 5 ++ 3 files changed, 358 insertions(+), 4 deletions(-) create mode 100644 datos/wines_2012.csv diff --git a/datos/wines_2012.csv b/datos/wines_2012.csv new file mode 100644 index 0000000..835503e --- /dev/null +++ b/datos/wines_2012.csv @@ -0,0 +1,181 @@ +judge,flight,wine,score,wine.amer,judge.amer +Jean-M Cardebat,white,A1,10,1,0 +Jean-M Cardebat,white,B1,13,1,0 +Jean-M Cardebat,white,C1,14,0,0 +Jean-M Cardebat,white,D1,15,0,0 +Jean-M Cardebat,white,E1,8,1,0 +Jean-M Cardebat,white,F1,13,1,0 +Jean-M Cardebat,white,G1,15,1,0 +Jean-M Cardebat,white,H1,11,0,0 +Jean-M Cardebat,white,I1,9,1,0 +Jean-M Cardebat,white,J1,12,0,0 +Tyler Colman,white,A1,16,1,1 +Tyler Colman,white,B1,14,1,1 +Tyler Colman,white,C1,14,0,1 +Tyler Colman,white,D1,16,0,1 +Tyler Colman,white,E1,12,1,1 +Tyler Colman,white,F1,11,1,1 +Tyler Colman,white,G1,11,1,1 +Tyler Colman,white,H1,14,0,1 +Tyler Colman,white,I1,11,1,1 +Tyler Colman,white,J1,14,0,1 +John Foy,white,A1,16,1,1 +John Foy,white,B1,17,1,1 +John Foy,white,C1,16,0,1 +John Foy,white,D1,15,0,1 +John Foy,white,E1,14.5,1,1 +John Foy,white,F1,14.5,1,1 +John Foy,white,G1,16,1,1 +John Foy,white,H1,17,0,1 +John Foy,white,I1,15,1,1 +John Foy,white,J1,17.5,0,1 +Olivier Gergaud,white,A1,14,1,0 +Olivier Gergaud,white,B1,19,1,0 +Olivier Gergaud,white,C1,12,0,0 +Olivier Gergaud,white,D1,10,0,0 +Olivier Gergaud,white,E1,19,1,0 +Olivier Gergaud,white,F1,18,1,0 +Olivier Gergaud,white,G1,17,1,0 +Olivier Gergaud,white,H1,16,0,0 +Olivier Gergaud,white,I1,18,1,0 +Olivier Gergaud,white,J1,14,0,0 +Robert Hodgson,white,A1,17,1,1 +Robert Hodgson,white,B1,11,1,1 +Robert Hodgson,white,C1,13,0,1 +Robert Hodgson,white,D1,14,0,1 +Robert Hodgson,white,E1,14,1,1 +Robert Hodgson,white,F1,10,1,1 +Robert Hodgson,white,G1,9,1,1 +Robert Hodgson,white,H1,9,0,1 +Robert Hodgson,white,I1,10,1,1 +Robert Hodgson,white,J1,10,0,1 +Linda Murphy,white,A1,15.5,1,1 +Linda Murphy,white,B1,15,1,1 +Linda Murphy,white,C1,17,0,1 +Linda Murphy,white,D1,18,0,1 +Linda Murphy,white,E1,16,1,1 +Linda Murphy,white,F1,17,1,1 +Linda Murphy,white,G1,15,1,1 +Linda Murphy,white,H1,14,0,1 +Linda Murphy,white,I1,16,1,1 +Linda Murphy,white,J1,17,0,1 +Daniele Meulder,white,A1,10,1,0 +Daniele Meulder,white,B1,15,1,0 +Daniele Meulder,white,C1,12,0,0 +Daniele Meulder,white,D1,12,0,0 +Daniele Meulder,white,E1,15,1,0 +Daniele Meulder,white,F1,14,1,0 +Daniele Meulder,white,G1,15,1,0 +Daniele Meulder,white,H1,12,0,0 +Daniele Meulder,white,I1,15,1,0 +Daniele Meulder,white,J1,12,0,0 +Jamal Rayyis,white,A1,16,1,0 +Jamal Rayyis,white,B1,15,1,0 +Jamal Rayyis,white,C1,14.5,0,0 +Jamal Rayyis,white,D1,17.5,0,0 +Jamal Rayyis,white,E1,16.5,1,0 +Jamal Rayyis,white,F1,14,1,0 +Jamal Rayyis,white,G1,12,1,0 +Jamal Rayyis,white,H1,15,0,0 +Jamal Rayyis,white,I1,13,1,0 +Jamal Rayyis,white,J1,12,0,0 +Francis Schott,white,A1,17,1,1 +Francis Schott,white,B1,16,1,1 +Francis Schott,white,C1,12,0,1 +Francis Schott,white,D1,18,0,1 +Francis Schott,white,E1,15,1,1 +Francis Schott,white,F1,16,1,1 +Francis Schott,white,G1,15,1,1 +Francis Schott,white,H1,14,0,1 +Francis Schott,white,I1,17,1,1 +Francis Schott,white,J1,15,0,1 +Jean-M Cardebat,red,A2,15,0,0 +Jean-M Cardebat,red,B2,11,0,0 +Jean-M Cardebat,red,C2,12,1,0 +Jean-M Cardebat,red,D2,16,1,0 +Jean-M Cardebat,red,E2,14,1,0 +Jean-M Cardebat,red,F2,11,1,0 +Jean-M Cardebat,red,G2,14.5,0,0 +Jean-M Cardebat,red,H2,13,1,0 +Jean-M Cardebat,red,I2,10,1,0 +Jean-M Cardebat,red,J2,14.5,0,0 +Tyler Colman,red,A2,14,0,1 +Tyler Colman,red,B2,11,0,1 +Tyler Colman,red,C2,16,1,1 +Tyler Colman,red,D2,12,1,1 +Tyler Colman,red,E2,14,1,1 +Tyler Colman,red,F2,13,1,1 +Tyler Colman,red,G2,14,0,1 +Tyler Colman,red,H2,12,1,1 +Tyler Colman,red,I2,13,1,1 +Tyler Colman,red,J2,11,0,1 +John Foy,red,A2,17.5,0,1 +John Foy,red,B2,19,0,1 +John Foy,red,C2,18,1,1 +John Foy,red,D2,18,1,1 +John Foy,red,E2,15,1,1 +John Foy,red,F2,16,1,1 +John Foy,red,G2,18,0,1 +John Foy,red,H2,18,1,1 +John Foy,red,I2,17,1,1 +John Foy,red,J2,17.5,0,1 +Olivier Gergaud,red,A2,10,0,0 +Olivier Gergaud,red,B2,17,0,0 +Olivier Gergaud,red,C2,9,1,0 +Olivier Gergaud,red,D2,14,1,0 +Olivier Gergaud,red,E2,19,1,0 +Olivier Gergaud,red,F2,12,1,0 +Olivier Gergaud,red,G2,15,0,0 +Olivier Gergaud,red,H2,10,1,0 +Olivier Gergaud,red,I2,11,1,0 +Olivier Gergaud,red,J2,18,0,0 +Robert Hodgson,red,A2,13,0,1 +Robert Hodgson,red,B2,17,0,1 +Robert Hodgson,red,C2,13,1,1 +Robert Hodgson,red,D2,16,1,1 +Robert Hodgson,red,E2,12,1,1 +Robert Hodgson,red,F2,15,1,1 +Robert Hodgson,red,G2,10,0,1 +Robert Hodgson,red,H2,12,1,1 +Robert Hodgson,red,I2,8,1,1 +Robert Hodgson,red,J2,11,0,1 +Linda Murphy,red,A2,13,0,1 +Linda Murphy,red,B2,14,0,1 +Linda Murphy,red,C2,17,1,1 +Linda Murphy,red,D2,16,1,1 +Linda Murphy,red,E2,15,1,1 +Linda Murphy,red,F2,17,1,1 +Linda Murphy,red,G2,14,0,1 +Linda Murphy,red,H2,15.5,1,1 +Linda Murphy,red,I2,13,1,1 +Linda Murphy,red,J2,18,0,1 +Daniele Meulder,red,A2,14,0,0 +Daniele Meulder,red,B2,16,0,0 +Daniele Meulder,red,C2,11,1,0 +Daniele Meulder,red,D2,16,1,0 +Daniele Meulder,red,E2,14,1,0 +Daniele Meulder,red,F2,15,1,0 +Daniele Meulder,red,G2,13,0,0 +Daniele Meulder,red,H2,11,1,0 +Daniele Meulder,red,I2,10,1,0 +Daniele Meulder,red,J2,15,0,0 +Jamal Rayyis,red,A2,15,0,0 +Jamal Rayyis,red,B2,19.5,0,0 +Jamal Rayyis,red,C2,14,1,0 +Jamal Rayyis,red,D2,12,1,0 +Jamal Rayyis,red,E2,13,1,0 +Jamal Rayyis,red,F2,16,1,0 +Jamal Rayyis,red,G2,14.5,0,0 +Jamal Rayyis,red,H2,15,1,0 +Jamal Rayyis,red,I2,16,1,0 +Jamal Rayyis,red,J2,16,0,0 +Francis Schott,red,A2,19,0,1 +Francis Schott,red,B2,18,0,1 +Francis Schott,red,C2,8,1,1 +Francis Schott,red,D2,15,1,1 +Francis Schott,red,E2,15,1,1 +Francis Schott,red,F2,12,1,1 +Francis Schott,red,G2,15,0,1 +Francis Schott,red,H2,16,1,1 +Francis Schott,red,I2,7,1,1 +Francis Schott,red,J2,17,0,1 diff --git a/notas/08-mcmc.qmd b/notas/08-mcmc.qmd index 7b90206..0ef3a47 100644 --- a/notas/08-mcmc.qmd +++ b/notas/08-mcmc.qmd @@ -621,7 +621,8 @@ y $H(\rho,\theta)$ es el Hamiltoniano en el paso anterior. **Observaciones**: 1. Un caso posible obtengamos desbordes o casi desbordes numéricos -del momento o la posición. Esto indica +del momento o la posición (el Hamiltoniano en el punto inicial es órdenes +de magnitud diferente que el inicial, ver [el manual de Stan](https://mc-stan.org/docs/reference-manual/mcmc.html#divergent-transitions) ). Esto indica problemas graves con el algoritmo de integración, y en general marcamos estas iteraciones como **divergentes**. Estas fallas pueden producir, como veremos, exploración insuficiente de la distribución objetivo. @@ -797,7 +798,6 @@ system.time(metropolis_2 <- metropolis_mc(1000, c(1,2), log_p, 1, 1)) #| eval: false #| fig-width: 10 #| fig-height: 5 -library(gganimate) sims_hmc <- hmc_1$sims |> mutate(n_sim = row_number()) |> mutate(algoritmo = "hmc") sims_metropolis_1 <- metropolis_1 |> @@ -831,14 +831,182 @@ la documentación de Stan. ## Diagnósticos de convergencia -Usualmente no es posible demostrar rigurosamente que un algoritmo MCMC es una -buena aproximación de la distribución posterior de interés- +Aunque casi nunca es posible demostrar rigurosamente que las simulaciones +de un algoritmo MCMC dan buena aproximación de la distribución posterior de interés, +especialmente con HMC y NUTS, tenemos muchos diagnósticos que fallan cuando +existen problemas serios. + +En primer lugar, será útil correr distintas cadenas con valores iniciales aleatorios diferentes, analizamos cada una y las comparamos entre sí. Recordamos que cada una +de estas cadenas tiene como distribución estacionaria límite la distribución posterior. +Diagnósticos que indican que las cadenas se comportan de manera muy distinta, explorando +distintas regiones del espacio de parámetros, o que +no han convergido porque exploran lentamente el espacio de parámetros, son señales +de problemas. + +Los diagnósticos más comunes son: + +1. Traza de cadenas +2. Medida R-hat de convergencia: mide la variabilidad entre cadenas y dentro de cadenas. +3. Número de muestras efectivas (ESS) y autocorrelación. +4. Transiciones divergentes. + + +### Modelos con variables latentes {-} + +Veremos el ejemplo de calificación de vinos de distintos países +de @rethinking, sus diagnósticos, +y aprovecharemos para introducir variables no observadas o latentes +para enriquecer nuestras herramientas de modelación. +Nuestra pregunta general es si el país de origen de los vinos influye en +su calidad. Los datos que tenemos son calificaciones de vinos de distintos +países por distintos jueces. La *calidad* del vino no la observamos directamente, +sino que es causa de las calificaciones que recibe. Para construir nuestro diagrama, +las consideraciones básicas son: +1. El origen del vino es una causa del calidad del vino (es nuestra cantidad a estimar). +2. Los jueces tienen distintas maneras de calificar, de manera que son causa de variación +en las calificaciones (hay jueces más duros, otros más barcos, etc.) No observamos directamente que tan "duro" es cada juez. +3. Los jueces califican vinos de distintos países de manera ciega. Sin embargo +es posible que reconozcan el país de origen por las características de los vinos, de +manera que puede existir un efecto directo de Origen en Calificación (no pasa por Calidad). +4. Es posible que Jueces de distintos países tienen distintos estándares de calificación. +```{r} +#| code-fold: true +grViz(" +digraph { + graph [ranksep = 0.2, rankdir = LR] + node [shape=plaintext] + Origen + Score + Origen_Juez + node [shape = circle] + Q + J + edge [minlen = 3] + Origen -> Q + Origen -> Score + Q -> Score + J -> Score + Origen_Juez -> J +} +") + +``` + +Y vemos, por nuestro análisis del DAG, que podemos identificar el efecto de Origen sobre +Calidad sin necesidad de estratificar por ninguna variable (no hay puertas traseras). +Sin embaergo, podemos estratificar por Juez para obtener más precisión (ver sección anterior de buenos y malos controles). +#### Primera iteración: modelo simple +Comenzamos con un modelo simple, y lo iremos construyendo para obtener la mejor estimación +posible de la influencia del país de origen en la calidad del vino. Nuestro primer modelo +consideramos que la calificación de cada vino depende de su calidad, y modelamos con una normal: + +$$S_i \sim \text{Normal}(\mu_i, \sigma)$$ +donde $$\mu_i = Q_{vino(i)}$$. Nuestra medida de calidad tiene escala arbitaria. Como +usaremos la calificación estandarizada, podemos poner +$$Q_j \sim \text{Normal}(0, 1).$$ +finalmente, ponemos una inicial para $\sigma$, por ejemplo $\sigma \sim \text{Exponential}(1)$ (puedes experimentar con una normal truncada también) + +```{r} +library(cmdstanr) +mod_vinos_1 <- cmdstan_model("./src/vinos-1.stan") +print(mod_vinos_1) +``` + +```{r} +# Wines 2022 de Statistical Rethinking +wines_2012 <- read_csv("../datos/wines_2012.csv") +glimpse(wines_2012) +wines_2012 <- wines_2012 |> + mutate(juez_num = as.numeric(factor(judge)), + vino_num = as.numeric(factor(wine))) |> + mutate(score_est = (score - mean(score))/sd(score)) +``` + + + +```{r} +n_jueces <- length(unique(wines_2012$juez_num)) +n_vinos <- length(unique(wines_2012$vino_num)) +c("num_vinos" = n_jueces, "num_jueces" = n_vinos, "num_datos" = nrow(wines_2012)) +``` +```{r} +datos_lst <- list( + N = nrow(wines_2012), + n_vinos = n_vinos, + n_jueces = n_jueces, + S = wines_2012$score_est, + vino = wines_2012$vino_num, + juez = wines_2012$juez_num +) +ajuste_vinos_1 <- mod_vinos_1$sample( + data = datos_lst, + chains = 4, + parallel_chains = 4, + iter_warmup = 1000, + iter_sampling = 2000, + refresh = 1000, + step_size = 0.1, +) + +``` + +Vemos que hay variabilidad en los vinos: + +```{r} +ajuste_vinos_1$summary(c("Q", "sigma")) |> + select(variable, mean, sd, q5, q95, rhat, ess_bulk, ess_tail) |> + filter(variable != "lp__") |> kable() +``` + + +Para hacer diagnósticos, podemos comenzar con las trazas de una cadena para todas las estimaciones de calidad de vino: + +```{r} +library(bayesplot) +mcmc_trace(ajuste_vinos_1$draws("Q", format = "df") |> filter(.chain == 1)) +``` +::: callout-tip +La traza de una cadena es la gráfica de las simulaciones de cada parámetro. Generalmente +buscamos que: no tenga tendencia, que no se quede "atorada" en algunos valores, y +que no muestre oscilaciones de baja frecuencia (la cadena "vaga" por los valores que explora). +::: + +Si incluímos todas las cadenas, nos fijemos en que todas ellas exploren +regiones similares del espacio de parámetros: + +```{r} +color_scheme_set("viridis") +mcmc_trace(ajuste_vinos_1$draws("Q", format = "df")) +``` +Lo que no queremos ver es lo siguiente, por ejemplo: + +```{r} +ajuste_vinos_malo <- mod_vinos_1$sample( + data = datos_lst, + chains = 4, + parallel_chains = 4, + iter_warmup = 5, + iter_sampling = 100, + refresh = 1000, + step_size =1 , + seed = 123 +) +``` + +```{r} +color_scheme_set("viridisA") +mcmc_trace(ajuste_vinos_malo$draws("Q", format = "df")) +``` +Hay varios problemas graves: +- Algunas cadenas parecen "atoradas" en ciertos valores +- Algunas cadenas parecen caminatas aleatorias (oscilaciones de baja frecuencia) +- Las cadenas no exploran de manera similar el espacio de parámetros diff --git a/renv.lock b/renv.lock index d632f27..0dd42f8 100644 --- a/renv.lock +++ b/renv.lock @@ -171,6 +171,11 @@ ], "Hash": "543776ae6848fde2f48ff3816d0628bc" }, + "bayesplot": { + "Package": "bayesplot", + "Version": "1.11.1", + "Source": "Repository" + }, "bit": { "Package": "bit", "Version": "4.0.5",