Es decir, el tamaño de los cálculos es una causa común de tratamiento (T) y resultado (M). Veremos más adelante que la decisión de condicionar a el tipo de cálculos proviene de un análisis relativamente simple de este diagrama causal, independientemente de los métodos que usemos para estimar las proporciones de interés (en este ejemplo, examinar las tablas cruzadas es equivalente a hacer estimaciones de máxima verosimlitud).
@@ -770,8 +776,8 @@Nótese que el análisis más apropiado no está en los datos: en ambos casos la tabla de datos es exactamente la misma. Los supuestos acerca del proceso que genera los datos sin embargo nos lleva a respuestas opuestas.
@@ -1226,6 +1232,24 @@Donde vemos ahora que el estado real de cada persona de la prueba es desconocido, aunque el resultado de la prueba depende de ese estado, y la cantidad de positivos que observamos es ahora \(N_{obs}\), que depende también de la sensibilidad y especificidad de la prueba.
@@ -551,8 +557,8 @@Usando argumentos como los del modelo original, las distribuciones de esp y sens son beta y podemos incorporarlas en la simulación de la posterior. Nuestra nueva función para simular el proceso generativo es:
@@ -646,7 +652,7 @@sims <- ajuste$draws(c("theta", "sens", "esp"), format = "df")
resumen <- ajuste$summary(c("theta"))Que también podríamos simplificar (suponiendo la \(N\) fija y conocida, pues \(N_+\) y \(M\) dan \(N_{-}\)) como:
@@ -361,8 +367,8 @@Y ahora construimos el modelo generativo. Supondremos que la muestra de \(N\) personas se toma de manera aleatoria de la población (una población grande, así que podemos ignorar el efecto de muestreo). Supondremos provisionalmente, además, que la prueba es perfecta, es decir, no hay falsos positivos o negativos.
@@ -1684,6 +1690,24 @@Omitiendo del diagrama las variables no observadas que también son causas únicamente de \(S\) y \(W, H\):
@@ -954,8 +960,8 @@Si queremos saber cómo influye el sexo en el peso, este diagrama indica que hay dos tipos de preguntas que podemos hacer:
@@ -1445,8 +1451,8 @@En este caso, el modelo causal es como sigue: conocemos la distancia \(D\) al hoyo en cada tiro. El éxito (\(Y=1\)) o fracaso (\(Y=0\)) depende de la distancia, junto con la velocidad a la que sale la pelota (muy alto o muy bajo puede dar un tiro fallido), y el ángulo \(\theta\) de salida. Adicionalmente, hay otros factors \(U\) que pueden afectar la probabilidad de éxito. Nótese que no escribiríamos, por ejemplo \(Y \leftarrow D\), porque la distancia no cambia causalmente con el resultado del tiro, aunque es cierto que si intervenimos en la distancia, esperaríamos obtener tasas de éxito diferentes. Igualmente, es necesario poner una flecha de \(V\) a \(D\) y \(V\) a \(Y\).
@@ -2383,6 +2389,24 @@Nótese que no describimos exactamente cómo son las funciones que relacionan las variables, sino más bien qué variables son causas directas de qué otras. Por ejemplo, aunque en nuestro ejemplo de arriba \(Y\) puede estar correlacionado con \(Z\), no hay una causa directa a \(Y\), porque cambios en \(Z\) afectan a \(X\), y es el cambio en \(X\) que es causa directa de \(Y\).
@@ -401,8 +407,8 @@En este ejemplos no podemos saber \(U1\) y \(U2\), y no nos interesa modelar la física de monedas, manera de lanzarlas, etc. En este ejemplo también no consideraremos qué hace que un día sea soleado o lluvioso (no nos interesa modelar el clima). En este momento, en teoría tenemos ecuaciones determinísticas para todas las variables, y si conocemos todas las variables exógenas \(U1,U2,U3,U4\) podríamos determinar exactamente lo que va a suceder con la ganancia, por ejemplo, o cualquier otra variable del sistema.
@@ -433,8 +439,8 @@# A tibble: 5 × 5
x d s1 s2 g
<dbl> <chr> <int> <int> <int>
-1 0.234 soleado 1 0 1
-2 0.281 soleado 3 3 3
-3 0.475 soleado 4 3 4
-4 0.944 lluvioso 5 5 10
-5 0.917 soleado 4 5 4En este caso,
@@ -627,8 +633,8 @@cor(sims_confusor |> select(x,y)) |> round(3) x y
-x 1.00 -0.43
-y -0.43 1.00Sin embargo, si condicionamos a \(Z\), que puede tomar los valores 0 o 1, vemos que \(X\) y \(Y\) son independientes, o dicho de otra manera, la condicional de \(Y\) dada \(Z\) y \(X\) sólo depende de \(Z\):
@@ -653,15 +659,15 @@cor(sims_confusor |> filter(z == 1) |> select(x,y)) |> round(3) x y
-x 1.000 -0.009
-y -0.009 1.000 x y
+x 1.000 0.008
+y 0.008 1.000cor(sims_confusor |> filter(z == 0) |> select(x,y)) |> round(3) x y
-x 1.000 -0.006
-y -0.006 1.000 x y
+x 1 0
+y 0 1Un ejemplo con variables continuas podría ser como sigue:
@@ -761,8 +767,8 @@Por la discusión de arriba, es claro que es necesario considerar la edad al casarse si queremos estimar el efecto de tasa de matrimonio en la tasa de divorcio. Es posible que la correlación entre estas dos tasas puede ser explicada solamente por la edad al casarse, y que en realidad al flecha \(M\to D\) sea muy débil o inexistente.
@@ -913,8 +919,8 @@Es decir, borramos todas las flechas que caen en \(M\) (pues la estamos interveniendo al valor que queramos), y luego simulando \(D\).
@@ -930,7 +936,7 @@En este caso,
@@ -1050,15 +1056,15 @@cor(sims_mediador |> filter(z == 1) |> select(x,y)) |> round(3) x y
-x 1.000 0.004
-y 0.004 1.000 x y
+x 1.000 -0.002
+y -0.002 1.000cor(sims_mediador |> filter(z == 0) |> select(x,y)) |> round(3) x y
-x 1.000 -0.002
-y -0.002 1.000Podemos también hacer un ejemplo continuo:
@@ -1151,8 +1157,8 @@Como el NSE es del hogar (una medida general de estatus social), se consideró en principio como una variable pre-tratamiento a la inteligencia de los niños por la que tradicionalmente se controlaba. Burks notó que hacer esto tenía no era apropiado, pues tiene como consecuencia cortar parte del efecto total de la inteligencia sobre el la inteligencia de los hijos. En otras palabras: la inteligencia de los padres hace más probable mejor NSE, y mejor NSE presenta mejores condiciones de desarrollo para sus hijos. Estatificar por esta variable bloquea este efecto.
@@ -1179,8 +1185,8 @@cor(sims_colisionador |> select(x,y)) x y
-x 1.000000000 0.008200524
-y 0.008200524 1.000000000Sin embargo, si condicionamos a \(Z\), que puede tomar los valores 0 o 1:
@@ -1253,8 +1259,8 @@cor(sims_colisionador |> filter(z == 0) |> select(x,y)) |> round(3) x y
-x 1.000 -0.271
-y -0.271 1.000print("Dado Z = 1")cor(sims_colisionador |> filter(z == 1) |> select(x,y)) |> round(3) x y
-x 1.000 0.377
-y 0.377 1.000Otro ejemplo con variables continuas:
@@ -1343,8 +1349,8 @@En este caso,
@@ -1370,8 +1376,8 @@[1] -0.003243749[1] 0.0003440308Sin embargo,
@@ -1396,16 +1402,16 @@cor(sims_colisionador |> filter(a ==0) |> select(x,y)) x y
-x 1.0000000 -0.2822012
-y -0.2822012 1.0000000cor(sims_colisionador |> filter(a ==1) |> select(x,y)) x y
-x 1.0000000 0.1097786
-y 0.1097786 1.0000000Consideremos la relación entre Z y Y. Primero vemos que hay dos caminos entre \(Z\) y \(Y\), que son \(p_1:X\gets V \to S\) y \(p_2: Z\to W \gets X \to Y\)
@@ -1687,8 +1693,8 @@Vimos que para calcular el efecto directo de \(F\) sobre \(W\), por ejemplo, es necesario bloquear el camino que pasa por \(G\) (estratificar por este nodo). Para el efecto total no es necesario condicionar a ningún otro nodo.
@@ -1816,8 +1822,8 @@En este caso:
@@ -1921,6 +1927,24 @@Nos interesa estimar el efecto causal de \(X\) sobre \(Y\). Sucede que en muchas ocasiones existen variables como \(U\) que son causas comunes de \(X\) y \(Y\). Como vimos, esto implica que no podemos simplemente ver la correlación entre \(X\) y \(Y\) para entender el efecto de \(X\) sobre \(Y\), pues una causa común de variación conjunta entre estas dos variables. Esta variable \(U\) puede ser observada o no.
@@ -350,8 +356,8 @@Nótese que:
@@ -409,8 +415,8 @@Gráfica de datos observacionales
Gráfica con intervención en A
Ahora queremos calcular \(p(a|do(d)) = p_m(a|d)\) en función de los datos. Siguiendo el mismo argumento que en el ejemplo anterior, sabemos que tenemos que estratificar o condicionar a \(T\) para poder usar nuestro proceso generador de observaciones, y obtenemos:
@@ -822,8 +828,8 @@Observamos que no podemos directamente usar la fórmula de ajuste pues NSE no es una variable observada.
@@ -852,8 +858,8 @@En este caso, todavía no podemos aplicar la fórmula original de ajuste pues no conocemos \(NSE\). Sin embargo, podemos bloquear los caminos no causales estratificando por Peso, y entonces podemos usar el criterio de puerta trasera para identificar el efecto del tratamiento, aún cuando no tengamos NSE.
@@ -1052,7 +1058,7 @@sims <- ajuste$draws( format = "df")
@@ -1173,7 +1179,7 @@ Ejemplo
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
-Chain 2 finished in 11.0 seconds.
+Chain 2 finished in 10.9 seconds.
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
@@ -1277,8 +1283,8 @@ ")
En este caso, el efecto de fumar (\(F\)) sobre cáncer (\(C\)) no es identificable pues no podemos condicionar a la variable de Genotipo (\(U\)). Supongamos que tenemos una medida adicional, que es la cantidad de depósitos de alquitrán den los pulmones de los pacientes. Este es es afectado por \(F\), y a su vez, el alquitrán incrementa la probabilidad de cáncer:
@@ -1304,17 +1310,17 @@La idea es primero estimar el efecto de \(F\) sobre \(A\), y después estimar el efecto de \(A\) sobre \(C\). La “composición” de estos dos efectos, dado el diagrama, debe darnos el estimador correcto. Primero consideramos el efecto de \(F\) sobre \(A\), y tenemos que (regla 2)
\[p(a|do(f)) = p(a|f),\] La igualdad se debe a que una vez que condicionamos a \(F\) no hay puertas traseras entre \(F\) y \(A\) (pues no condicionamos a \(C\)). Esta dependencia causal la podemos entonces estimar de los datos.
-El efecto de \(A\) sobre \(C\) también es identificable, pues el camino no causal se bloquea cuando condicionamos a \(A\), de forma que por la fórmula de ajuste:
+El efecto de \(A\) sobre \(C\) también es identificable, pues el camino no causal se bloquea cuando condicionamos a \(F\), de forma que por la fórmula de ajuste:
\[p(c|do(a)) = \int p(c|a, f') p(f')\, df'\]
Ahora encadenamos estas dos ecuaciones:
\[p(c|do(f)) = \int p(c|do(a))p(a|f)\,da\]
-que equivale en simulación a: dado un valor de \(F\), simulamos \(A=a\) con nuestro modelo ajustado con datos naturales. Ahora intervenimos \(A\) con el valor a que obtuvimos y simulamos \(C\). Sin embargo, para hacer este último paso con datos naturales, necesitamos usar el criterio de puerta trasera como explicamos arriba: simulamos entonces \(f´\) de \(p(f)\), y después simulamos \(C\) en función de \(a\) y \(f´\) (con una distribución construida a partir de datos).
+que equivale en simulación a: dado un valor de \(F\), simulamos \(A\) con nuestro modelo ajustado con datos naturales. Ahora intervenimos \(A\) con el valor \(a\) que obtuvimos y simulamos \(C\). Sin embargo, para hacer este último paso con datos naturales, necesitamos usar el criterio de puerta trasera como explicamos arriba: simulamos entonces \(f´\) de \(p(f)\), y después simulamos \(C\) en función de \(a\) y \(f´\) (con una distribución construida a partir de datos).
Requerimos en este caso construir y estimar la condicional \(p(c|a, f)\) basado en los datos.
En fórmula, en general, se escribe como:
+ +No hay ninguna variable confusora, y una estrategia de estimación es comparar \(PF\) entre los grupos.
@@ -492,8 +499,8 @@Hemos añadido un nodo implícito (otros factores que afectan \(Y\) y no tienen relación con otras variables del sistema) para explicar qué es lo que pasa cuando condicionamos a \(Z\): como \(Z\) es un descendiente del colisionador en \(Y\), se activa una ruta no causal entre \(U_y\) y \(T\), y estas dos cantidades aparecen como correlacionadas (es una correlación no causal). Esto en consecuencia modifica la correlación entre \(T\) y \(Y\).
@@ -622,8 +629,8 @@En la gráfica de arriba, \(T\) indica si la madre es fumadora o no, y \(Y\) la mortalidad. \(Z\) si el bebé nació con bajo peso o no.
@@ -681,8 +688,8 @@En este caso, condicionar a \(Z\) no sesga nuestras estimaciones, pues no activamos ninguna ruta no causal. La dificultad es que típicamente disminuye la precisión de la estimación (usamos un modelo más grande donde no es necesario):
@@ -728,8 +735,8 @@En este caso, tenemos una variable confusora \(U\) que no nos permite estimar sin sesgo el efecto de \(T\) sobre \(Y\). Sin embargo, si condicionamos a \(Z\), la situación puede emperorar (amplificación de sesgo), pues dentro de cada nivel de \(Z\) hay menos variación de \(X\), y eso implica que la covarianza entre \(X\) y \(Y\), en cada nivel de \(Z\), se debe más a la variable confusora.
@@ -855,6 +862,24 @@library(tidyverse)
+library(kableExtra)
+library(DiagrammeR)
+ggplot2::theme_set(ggplot2::theme_light())
+inv_logit <- \(x) 1 / (1 + exp(-x)) En esta sección exlicaremos brevemente cómo funcionan paquetes como Stan para producir simulaciones de una posteriores complicadas en dimensión alta.
+En primer lugar, recordemos que si queremos calcular la posterior de un modelo (generalmente para calcular después resúmenes que involucran integrales de esta posterior) tenemos los siguiente enfoques:
+Stan utiliza 3, y hay variedad de algoritmos MCMC. Ya discutimos que 1, la aproximación analítica, es en general imposible (a menos fuera de ciertos modelos restringidos). La aproximación 2 excesivamente intensiva, al grado que sólo para modelos muy chicos y con pocos parámetros es posible utilizarla. Existen otros métodos también como aproximaciones cuadráticas que en ciertos casos funcionan, pero son limitados en su aplicación.
+La idea de simulación de variables aleatorias es ahora fundamental en muchas áreas científicas, incluyendo la estadística, y los métodos que la utilizan se llaman métodos de Monte Carlo. Por ejemplo, considera el método bootstrap, pruebas de permutaciones, validación cruzada, y en general simulación para cálculo de resúmenes que son difíciles de calcular directamente (por ejemplo, la mediana de una distribución Gamma, ver Median approximations and bounds aquí).
+Uno de los primeros algoritmos MCMC fue el de Metropolis-Hastings, que veremos primero en algunos ejemplos. Veremos también por qué ahora tenemos mejores opciones que MH para estimar posteriores de nuestros modelos.
+Supongamos que queremos simular de una variable aleatoria \(X\) con distribución discreta sobre los valores \(1,2\ldots, k\) con probabilidades \(p(1),p(2),\ldots,p(k)\). Este problema puede resolverse fácilmente de varias maneras, pero utilizaremos un método de Monte Carlo tipo Metropolis. Supongamos que podemos simular de una variable aleatoria \(U\) que es uniforme en \(1,2,\ldots, k\) (con probabilidades iguales a 1/k).
+Lo que podemos hacer es lo que sigue, para \(i=1,\ldots, M\):
+Para cada \(i\),
+El resultado es una sucesión de valores \(x_1,x_2,\ldots, x_M\). Es posible demostrar que la distribución de estas \(x_i\) converge a la distribución \(p(1),\ldots, p(k)\) si \(M\) es suficientemente grande.
+Este método se llama Metropolis-Hastings. Es un método de Monte Carlo, y como podemos ver, se trata de una cadena de Markov, pues la distribución cada siguiente lugar \(x_{i+1}\), condicionada al valor actual \(x_i\) no depende de valores anteriores de las \(x\).
+#set.seed(451123)
+# definimos estas p
+k <- 40
+p <- exp(-(1:k - k/3)^2 / 10)
+p <- p /sum(p)
+dist_obj <- tibble(x = 1:k, p = p)
+# simulamos
+M <- 200000
+x <- numeric(M)
+x[1] <- 20
+for(i in 1:M){
+ u <- sample(1:k, 1)
+ q <- p[u] / p[x[i]]
+ if(runif(1) < q){
+ x[i+1] <- u
+ } else {
+ x[i+1] <- x[i]
+ }
+}En rojo mostramos las probabilidades objetivo que queremos estimar, y en negro las estimadas con nuestro método de arriba.
+resultados_sim <- tibble(x = x) |>
+ mutate(n_sim = row_number())
+resumen_sim <- resultados_sim |> count(x) |>
+ right_join(tibble(x = 1:k, p = p)) |>
+ mutate(n = ifelse(is.na(n), 0, n)) |>
+ mutate(p_aprox = n / sum(n))Joining with `by = join_by(x)`ggplot(dist_obj, aes(x = x)) +
+ geom_point(aes(y = p)) +
+ geom_point(data = resumen_sim,
+ aes(y = p_aprox), color = "red", size = 3, alpha = 0.5) 
Como vemos, los valores de \(x_1,\ldots, x_M\) se distribuyen aproximadamente como la distribución \(p\) objetivo. Esta es una manera de simular valores de esta distribución discreta \(p\). Podemos ver cómo se ven las simulaciones sucesivas:
+ggplot(resultados_sim |> filter(n_sim < 400), aes(x = n_sim, y = x)) +
+ geom_line() + scale_y_continuous(breaks = 1:20)
El defecto que tiene este algoritmo es que puede ser relativamente lento, pues vemos que hay periodos largos donde se “atora” en valores de probabilidad relativamente alta. La razón es que en muchos pasos, estamos proponiendo “saltos al vacío” a lugares de probabilidad muy baja, que rara vez se aceptan.
+Podemos hacer más eficiente nuestro algoritmo si le permitimos explorar con mayor facilidad los posibles valores de \(x\). Esto se logra proponiendo saltos locales: si estamos en \(x_i\), entonces proponemos los valores \(x_i - 1\) y \(x_i + 1\) con la misma probabilidad 1/2 (excepto en los extremos donde sólo proponemos \(x_i,x_i+1\) o \(x_i-1,x_i\)).
+Proponemos entonces la suguiente modificación del paso 1 de propuesta:
+Para cada \(i\),
+Esto lo escribimos como sigue:
+#set.seed(4511)
+# simulamos
+x <- numeric(M)
+x[1] <- 20
+for(i in 1:M){
+ u <- sample(c(x[i] - 1, x[i] + 1), 1)
+ if(u == k+1) u <- k
+ if(u == 0) u <- 1
+ q <- p[u] / p[x[i]]
+ if(runif(1) < q){
+ x[i+1] <- u
+ } else {
+ x[i+1] <- x[i]
+ }
+}Obtenemos:
+resultados_sim_2 <- tibble(x = x) |>
+ mutate(n_sim = row_number())
+resumen_sim_2 <- resultados_sim_2 |> count(x) |>
+ mutate(p_aprox = n / sum(n))
+ggplot(dist_obj, aes(x = x)) +
+ geom_point(aes(y = p)) +
+ geom_point(data = resumen_sim_2,
+ aes(y = p_aprox), color = "red", size = 3, alpha = 0.5) 
Y podemos ver cómo evoluciona nuestra cadena de Markov:
+ggplot(resultados_sim_2 |> filter(n_sim < 400), aes(x = n_sim, y = x)) +
+ geom_line() + scale_y_continuous(breaks = 1:10)
¿Cómo se comparan estos dos métodos? Podemos ver por ejemplo cómo se comparan las distribuciones aproximadas hasta cierto número de iteraciones con la verdadera distribución objetivo:
+approx_sim <- map_df(seq(200, 30000, by = 200), function(n_sims){
+ resumen_1 <- resultados_sim |> filter(n_sim <= n_sims) |>
+ count(x) |>
+ mutate(p_aprox = n / sum(n)) |>
+ select(-n) |>
+ right_join(dist_obj, by = "x") |>
+ mutate(metodo = "MH-1") |>
+ mutate(n_sims = n_sims)
+ resumen_2 <- resultados_sim_2 |> filter(n_sim <= n_sims) |>
+ count(x) |>
+ mutate(p_aprox = n / sum(n)) |>
+ select(-n) |>
+ right_join(dist_obj, by = "x") |>
+ mutate(metodo = "MH-2") |>
+ mutate(n_sims = n_sims)
+ bind_rows(resumen_1, resumen_2) |>
+ mutate(p_aprox = ifelse(is.na(p_aprox), 0, p_aprox))
+})approx_sim |>
+ mutate(dif_abs = abs (p_aprox-p)) |>
+ group_by(metodo, n_sims) |>
+ summarise(dif_abs = sum(dif_abs)) |>
+ggplot(aes(n_sims, dif_abs, color = metodo)) +
+ geom_line() `summarise()` has grouped output by 'metodo'. You can override using the
+`.groups` argument.
En este caso, considera qué es lo que sucede en cada uno de estos casos:
+Este es el primer balance que existe en este algoritmo: tomar pasos grandes y balancear las probabilidades quizá rechazando muy frecuentemente (no es eficiente), o tomar pasos chicos y vagar más tiempo para visitar regiones de alta probabilidad, aunque con menos tasa de rechazo. Dependiendo de la distribución que queremos aproximar podemos inclinarnos más por una o por otra opción.
+¿Por qué funcionan estos algoritmos? Supongamos que en cada paso, se cumple que (balance detallado): \[{q(x|y)}p(y) = {q(y|x)}{p(x)}\] donde \(q(x|y)\) son las probabilidades de transición de nuestra cadena de Markov propuesta. Esta ecuación dice que si la probabilidad se distribuye como \(p(x)\), entonces al transicionar con \(q\) la probabilidad fluje de manera que se mantiene estática en \(p\), es decir \(p\) es una distribución estacionaria. Esto es fácil de demostrar pues \[\sum_{y} q(x|y)p(y) = \sum_{y} q(y|x)p(x) = p(x) \sum_{y} q(y|x) = p(x).\]
+Bajo otros supuestos adicionales de ergodicidad (aperiodicidad y tiempos de recurrencia finitos, es decir, las transiciones mezclan bien los estados), entonces podemos simular la cadena de Markov por un tiempo suficientemente largo y con esto obtener una muestra de la distribución objetivo \(p\).
+¿Cómo podemos diseñar entonces las \(q(x|y)\) correspondientes? Comenzamos considerando distribuciones propuesta \(q_0(x|y)\) que no satisfacen la ecuación de balance, y supondremos como en los ejemplos de arriba (verifícalo) que nuestras transiciones tienen probabilidades simétricas \(q_0(y|x) = q_0(x|y)\). Entonces, cuando \(p(y)/p(x) > 1\), queremos transicionar de \(x\) a \(y\) con más frecuencia que de \(y\) a \(x\). Comenzando en \(x\), si la propuesta de \(q_0\) es \(y\), podríamos poner entonces que el sistema transicione con probabilidad 1 a \(y\). Sin embargo, si empezamos en \(y\) y la propuesta es \(x\), el sistema sólo transiciona de \(y\) a \(x\) con probabilidad \(p(x)/p(y)\).
+De esta manera, obtenemos que bajo \(q(y|x)\), \(x\) transiciona a \(y\) con probabilidad \(\min\{1, p(y)/p(x)\}\). Entonces, el cociente \(\frac{q(y|x)}{q(x|y)}\) es igual a \(\frac{p(y)}{p(x)}\) si \(p(y)<p(x)\), y es igual a \(1/\frac{p(x)}{p(y)} = \frac{p(y)}{p(x)}\) si \(p(y)>p(x)\).
+En MCMC, buscamos un cadena de Markov que, en el largo plazo, visite cada posible parámetro proporcionalmente a la probabildad posterior de cada parámetro.
+Supongamos ahora que quisiéramos simular de una normal multivariada con media en c(2,3) y matriz de covarianza \(\Sigma\), que supondremos es tal que la desviación estándar de cada variable es 1 y la correlación es 0.8. La matriz \(\Sigma\) tiene 1 en la diagonal y 0.8 fuera de la diagonal.
+La distribución objetivo \(p\) está dada entonces (módulo una constante de proporcionalidad):
+Sigma <- matrix(c(1, 0.8, 0.8, 1), nrow = 2)
+Sigma_inv <- solve(Sigma)
+m <- c(2, 3)
+log_p <- function(z){
+ - 0.5 * (t(z-m) %*% Sigma_inv %*% (z-m))
+}Nótese que como Metropolis hace cocientes de probabilidades, sólo es necesario conocer la densidad objetivo módulo una constante de proporcionalidad.
+Una algoritmo de Metropolis podría ser el siguiente:
+# simulamos
+M <- 50000
+simular_normal <- function(M, m, Sigma, delta_x, delta_y){
+ z <- matrix(nrow = M, ncol = 2)
+ z[1, ] <- c(2.5, 3.2)
+ colnames(z) <- c("x", "y")
+ rechazo <- 0
+ for(i in 1:(M-1)){
+ x_prop <- rnorm(1, z[i, 1], delta_x)
+ y_prop <- rnorm(1, z[i, 2], delta_y)
+ z_prop <- c(x_prop, y_prop)
+ q <- exp(log_p(z_prop) - log_p(z[i, ]))
+ if(runif(1) < q){
+ z[i + 1, ] <- z_prop
+ } else {
+ rechazo <- rechazo + 1
+ z[i + 1, ] <- z[i, ]
+ }
+ }
+ print(rechazo / M)
+ z_tbl <- as_tibble(z) |>
+ mutate(n_sim = row_number())
+ z_tbl
+}
+z_tbl <- simular_normal(M, m, Sigma, 1, 1)[1] 0.60038Vemos que tenemos una tasa alta de rechazos. ¿Por qué? Veamos cómo se ven las simulaciones después de 1000 iteraciones:
+# estas las usamos para graficar
+sims_normal <- mvtnorm::rmvnorm(10000, mean = m, sigma = Sigma)
+colnames(sims_normal) <- c("x", "y")
+sims_normal <- as_tibble(sims_normal)
+graf_tbl <- map_df(seq(10, 500, 20), function(i){
+ z_tbl |> filter(n_sim <= i) |> mutate(num_sims = i)
+})
+ggplot(graf_tbl, aes(x, y)) +
+ stat_ellipse(data = sims_normal, aes(x, y),
+ level = c( 0.9), type = "norm", colour = "salmon") +
+ stat_ellipse(data = sims_normal, aes(x, y),
+ level = c( 0.5), type = "norm", colour = "salmon") +
+ geom_point(alpha = 0.1) +
+ facet_wrap(~num_sims) + theme_minimal()
Observaciones:
+Podemos entonces proponer saltos más chicos, por ejemplo:
+z_tbl <- simular_normal(M, m, Sigma, 0.2, 0.2)[1] 0.1568graf_tbl <- map_df(seq(10, 500, 20), function(i){
+ z_tbl |> filter(n_sim <= i) |> mutate(num_sims = i)
+})
+ggplot(graf_tbl, aes(x, y)) +
+ stat_ellipse(data = sims_normal, aes(x, y),
+ level = c( 0.9), type = "norm", colour = "salmon") +
+ stat_ellipse(data = sims_normal, aes(x, y),
+ level = c( 0.5), type = "norm", colour = "salmon") +
+ geom_point(alpha = 0.1) +
+ facet_wrap(~num_sims) + theme_minimal()
Observaciones:
+En el algoritmo de Metropolis Hastings hay una tensión natural en tamaño de salto y tasa de aceptación. Si el tamaño de los saltos es muy grande, la tasa de aceptación es baja y esto producen ineficiencias. Si el tamaño de los saltos es muy chico, la tasa de aceptación es más alta, pero esto también es ineficiente.
+Existen métodos que pueden superar este problema, como son muestreo de Gibbs y Monte Carlo Hamiltoniano. El primero no lo discutiremos, pues requiere poder simular fácilmente de cada parámetro dados los otros, y esto no siempre es posible. Veremos más el segundo, donde usaremos información del gradiente de la distribución objetivo para proponer exploración más eficiente.
+