Agregar principio de diagnósticos

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

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

renv.lock

@@ -171,6 +171,11 @@
       ],
       "Hash": "543776ae6848fde2f48ff3816d0628bc"
     },
+    "bayesplot": {
+      "Package": "bayesplot",
+      "Version": "1.11.1",
+      "Source": "Repository"
+    },
     "bit": {
       "Package": "bit",
       "Version": "4.0.5",