¿Has estudiado estadística bayesiana?

  1. Sí
  2. No

En las siguientes dos clases estudiaremos métodos bayesianos, con el fin de entender las diferencias entre el paradigma bayesiano y frecuentista comenzamos recordando algunos aspectos de cada uno. El punto de vista frecuentista se basa en los siguientes puntos (Wasserman, 2005):

  1. La probabilidad se refiere a un límite de frecuencias relativas, 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%.

Por su parte, el paradigma Bayesiano se basa en 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 copa de te el 1ro 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 infernecias como estimaciones puntuales y estimaciones de intervalos se pueden extraer de dicha distribución.

Probabilidad subjetiva

¿Qué tanta certeza tienes de que una moneda acuñada por la casa de moneda mexicana es justa? Si, en cambio, consideramos una moneda antigua y asimétrica, ¿creemos que es justa? En estos escenarios no estamos considerando la verdadera probabilidad, inherente a la moneda, lo que queremos medir es el grado en que creemos que cada probabilidad puede ocurrir.

Con el fin de especificar nuestras creencias debemos medir que tan verosimil pensamos que es cada posible resultado. Por ejemplo, puedes pensar que una mujer mexicana promedio mide 156 cm pero estar abierto a la posibilidad de que el promedio sea un poco mayor o menor. Es así que puedes describir tus creencias a través de una curva con forma de campana y centrada en 156. No olvidemos que estamos describiendo probabilidades y por tanto deben cumplir los axiomas de probabilidad. Es por esto que la curva debe conformar una distribuión de probabilidad.

Si \(p(\theta)\) representa el grado de nuestra creencia en los valores de \(\theta\), entonces la media de \(p(\theta)\) se puede pensar como un valor de \(\theta\) que representa nuestra creencia típica o central. Por su parte, la varianza de \(\theta\), que mide que tan dispersa esta la distribución, se puede pensar como la incertidumbre entre los posibles valores.

Regla de Bayes

La regla de Bayes tiene un papel fundamental en la estadística bayesiana, se utiliza para determinar la probabilidad de un modelo dado un conjunto de datos. Recordemos que lo que el modelo determina es la probabilidad de los datos condicional a valores particulares de los parámetros y a la estructura del modelo. Por otra parte usamos la regla de Bayes para ir de la probabilidad de los datos, dado el modelo, a la probabilidad del modelo, dados los datos.

Ejemplo: Lanzamientos de monedas

Comencemos recordando la regla de Bayes usando dos variables aleatorias discretas. Lanzamos una moneda 3 veces, sea X la variable aleatoria correspondiente al número de Águilas observadas y Y registra el número de cambios entre águilas y soles.

  • Escribe la distribución conjunta de las variables.

  • Considera la probabilidad de observar un cambio condicional a que observamos un águila y compara con la probabilidad de observar un águila condicional a que observamos un cambio.

Para entender probabilidad condicional podemos pensar en restringir nuestra atención a una única fila o columna de la tabla. Supongamos que alguien lanza una moneda 3 veces y nos informa que la secuencia contiene exactamente un cambio. Dada esta información podemos restringir nuestra atención a la fila corrspondiente a un solo cambio. Sabemos que ocurrió uno de los eventos de esa fila. Las probabilidades relativas de los eventos de esa fila no han cambiado pero sabemos que la probabilidad total debe sumar uno, por lo que simplemente normalizamos dividiendo entre \(p(C=1)\).

Regla de Bayes en modelos y datos

Ahora consideremos el caso en que las variables fila y columna corresponden a datos y parámetros del modelo respectivamente. Un modelo especifica la probabilidad de valores particulares dado la estructura del modelo y valores de los parámetros. Por ejemplo en un modelo de lanzamientos de monedas tenemos \(p(x = A|\theta)=\theta\) y \(p(x = S|\theta)= 1- \theta\).

De manera general, el modelo especifica: \[p(datos|\text{valores de parámetros y estructura del modelo})\]

y usamos la regla de Bayes para convertir la expresión anterior a lo que nos interesa de verdad, que es, que tanta certidumbre tenemos del modelo condicional a los datos:

\[p(\text{valores de parámetros y estructura del modelo} | datos)\]

Una vez que observamos los datos, usamos la regla de Bayes para determinar o actualizar nuestras creencias de posibles parámetros y modelos.

Antes de proceder hace falta definir algunos conceptos y notación: \[p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)}\]

\(p(\theta|x)\) se conoce como la distribución posterior, \(p(x|\theta)=\mathcal{L}(\theta)\) es la verosimilitud, \(p(\theta)\) es la distribución inicial o a priori y \(p(x)\) es la evidencia.

  • La distribución inicial \(p(\theta)\) cuantifica la información de \(\theta\) externa a \(x\), esto es, la distribución a priori resume nuestras creencias acerca del parámetro ajenas a los datos.

  • La distribución de verosimilitud cuantifica la información de \(\theta\) asociada a los datos.

  • La evidencia \(p(x)\) es la probabilidad de los datos de acuerdo al modelo, se determina sumando (o integrando) a través de todos los posibles valores de los parámetros, ponderado por la certeza en dichos parámetros.

Calculo de la distribución posterior

En la inferencia Bayesiana se requiere calcular el denominador de la fórmula de Bayes \(p(x)\), es común que esto requiera que se calcule una integral complicada; sin embargo hay algunas maneras de evitar esto,

  1. El camino tradicional consiste en usar funciones de verosimilitud con dsitribuciones iniciales conjugadas. Cuando una distribución inicial es conjugada de la verosimilitud resulta en una distribución posterior con la misma forma funcional que la distribución inicial.

  2. Otra alternativa es aproximar la integral numericamente. Cuando el espacio de parámetros es de dimensión chica, se puede cubrir con una cuadrícula de puntos y la integral se puede calcular sumando a través de dicha cuadrícula. Sin embargo cuando el espacio de parámetros aumenta de dimensión el número de puntos necesarios para la aproximación crece demasiado y hay que recurrir a otas técnicas.

  3. Se ha desarrollado una clase de métodos de simulación para poder calcular la distribución posterior, estos se conocen como cadenas de Markov via Monte Carlo (MCMC por sus siglas en inglés). El desarrollo de los métodos MCMC es lo que ha propiciado el desarrollo de la estadística bayesiana en años recientes.

Ejemplo: Bernoulli

Distribuciones conjugadas

Comenzaremos con el modelo Beta-Binomial. Recordemos que si X en un experimento con dos posibles resultados, X se distribuye Bernoulli y la función de probabilidad esta definida por:

\[p(x|\theta)=\theta^x(1-\theta)^{1-x}\]

recordemos también que podemos pensar en la fórmula anterior como una función de \(x\), donde \(x\) toma uno de dos valores 0 o 1. También podemos pensar que \(x\) esta fija (observada) y la expresión es una función de \(\theta\), en este caso, la ecuación especifica la probabilidad de un valor fijo \(x\) para un valor de \(\theta\) y la llamamos la verosimilitud de \(\theta\). En inferencia Bayesiana la función \(p(x|theta)\), usualmente se considera como función de \(\theta\).

Ahora, si lanzamos una moneda \(N\) veces tenemos un conjunto de datos \(\{x_1,...,x_N\}\), suponemos que los lanzamientos son independientes por lo que la probabilidad de observar el conjunto de \(N\) lanzamientos es el producto de las probabilidades para cada observación:

\[p(x_1,...,x_N|\theta) = \prod_{n=1}^N p(x_n|\theta)\] \[= \theta^z(1 - \theta)^{N-z}\]

donde \(z\) denota el número de éxitos (águilas).

Ahora, en principio para describir nuestras creencias iniciales podríamos usar cualquier función de densidad con soporte en \([0, 1]\), sin embargo, sería conveniente que el producto \(p(x|\theta)p(\theta)\) (el numerador de la fórmula de Bayes) resulte en una función con la misma forma que \(p(\theta)\). Cuando este es el caso, las creencias inicial y posterior se describen con la misma distribución. Esto es conveninte pues si obtenemos nueva información podemos actualizar nuestro conocimiento de manera inmediata, conservando la forma de las distribuciones.

Cuando las funciones \(p(x|\theta)\) y \(p(\theta)\) se combinan de tal manera que la distribución posterior pertenece a la misma familia (tiene la misma forma) que la distribución inicial, entonces decimos que \(p(\theta)\) es conjugada para \(p(x|\theta)\).


Vale la pena notar que la inicial es conjugada únicamente respecto a una función de verosimilitud particular.

Una distribución conjugada para \(p(x|\theta) = \theta^z(1 - \theta)^{N-z}\) es una Beta(a, b) \[p(\theta) = \frac {\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}\]

Para describir nuestro conocimiento inicial podemos explorar la media y desviación estándar de la distribución beta, la media es \[\bar{\theta} = a/(a+b)\] por lo que si \(a=b\) la media es 0.5 y conforme aumenta \(a\) en relación a \(b\) aumenta la media. La desviación estándar es

\[\sqrt{\bar{\theta}(1-\bar{\theta})/(a+b+1)}\]

Una manera de seleccionar los parámetros \(a\) y \(b\) es pensar en la proporción media de águilas (m) y el tamaño de muestra (n). Ahora, \(m=a/(a+b)\) y \(n = a+b\), obteniendo.

\[a=mn, b=(1-m)n\]

Otra manera es comenzar con la media y la desviación estándar. Al usar este enfoque debemos recordar que la desviación estándar debe tener sentido en el contexto de la densidad beta. En particular la desviación estándar típicamente es menor a 0.289 que corresponde a la desviación estándar de una uniforme. Entonces, para una densidad beta con media \(m\) y desviación estándar \(s\), los parámetros son: \[a=m\bigg(\frac{m(1-m)}{s^2}- 1\bigg), b=(1-m)\bigg(\frac{m(1-m)}{s^2}- 1\bigg)\]

Una vez que hemos determinado una inicial conveniente para la verosimilitud Bernoulli calculamos la posterior. Supongamos que observamos \(N\) lanzamientos de los cuales \(z\) son águilas, entonces podemos ver que la posterior es nuevamente una densidad Beta.

\[p(\theta|z)\propto \theta^{a+z-1}(1 -\theta)^{(N-z+b)-1}\]

Concluímos entonces que si la distribución inicial es Beta(a,b), la posterior es Beta(z+a,N-z+b).

Vale la pena explorar la relación entre la distribución inicial y posterior en las medias. La media incial es \[a/(a+b)\] y la media posterior es \[(z+a)/[(z+a) + (N-z+b)]=(z+a)/(N+a+b)\] podemos hacer algunas manipulaciones algebráicas para escribirla como:

\[\frac{z+a}{N+a+b}=\frac{z}{N}\frac{N}{N+a+b} + \frac{a}{a+b}\frac{a+b}{N+a+b}\]

es decir, podemos escribir la media posterior como un promedio ponderado entre la media inicial \(a/(a+b)\) y la proporción observada \(z/N\).

Ahora podemos pasar a la inferencia, comencemos con estimación de la proporción \(\theta\). La distribución posterior resume todo nuestro conocimiento del parámetro \(\theta\), en este caso podemos graficar la deistribución y extraer valores numéricos como la media.

grid.xy <- data.frame(x = seq(0, 1, .1), y = seq(0, 1, .1))
# N = 14, z = 11, a = 1, b = 1
ggplot(grid.xy, aes(x = x, y = y)) + 
  stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 1), 
    color = "orange") +  # inicial
  stat_function(fun = dbeta, args = list(shape1 = 12, shape2 = 4), 
    color = "green3") +  # verosimilitud
  stat_function(fun = dbeta, args = list(shape1 = 12, shape2 = 4), 
    color = "steelblue") # posterior

# N = 14, z = 11, a = 100, b = 100 
ggplot(grid.xy, aes(x = x, y = y)) + 
  stat_function(fun = dbeta, args = list(shape1 = 100, shape2 = 100), 
    color = "orange") + 
  stat_function(fun = dbeta, args = list(shape1 = 12, shape2 = 4), 
    color = "green3") + 
  stat_function(fun = dbeta, args = list(shape1 = 111, shape2 = 103), 
    color = "steelblue")

Una manera de resumir la distribución posterior es a través de intervalos de probabilidad, otro uso de los intervalos es establecer que valores del parámetro son creíbles.

Calcula un intervalo del 95% de probabilidad para cada una de las distribuciones posteriores del ejemplo.

Ahora pasemos a predicción, calculamos la probabilidad de \(y =1\):

\[p(y = 1) = \int p(y=1|\theta)p(\theta|z)d\theta\] \[=\int \theta p(\theta|z,N) d\theta\] \[=(z+a)/(N+a+b)\] Esto es, la probabilidad predictiva de águila es la media de la distribución posterior sobre \(\theta\).

Metrópolis

Hay ocasiones en las que los métodos de inicial conjugada y aproximación por cuadrícula no funcionan, hay casos en los que la distribución beta no describe nuestras creencias iniciales. Por su parte, la aproximación por cuadrícula no es factible cuando tenemos varios parámetros. Es por ello que surge la necesidad de utilizar métodos de Monte Carlo vía Cadenas de Markov (MCMC).

En Metropolis se asume que podemos calcular \(p(\theta)\) para un valor particular de \(\theta\) y el valor de la verosimilitud \(p(x|\theta)\) para cualquier \(x\), \(\theta\) dados. En realidad, el método únicamente requiere que se pueda calcular el producto de la inicial y la verosimilitud, y sólo hasta una constante de proporcionalidad. Lo que el método produce es una aproximación de la distribución posterior \(p(\theta|x)\) mediante una muestra de valores \(\theta\) obtenido de dicha distribución.

Caminata aleatoria. Con el fin de entender el algoritmo comenzaremos estudiando el concepto de caminata aleatoria. Supongamos que un vendedor de enciclopedias 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; 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_{prop}\) la población de la isla propuesta y \(P_{actual}\) la población de la isla actual. Entonces el vendedor cambia de isla con probabilidad \(p_{mover}=P_{prop}/P_{actual}\).

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.

islas <- data.frame(islas = 1:10, pob = 1:10)

caminaIsla <- 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 estas en los extremos y el volado indica salir
    return(i)
  }
  u2 <- runif(1)
  p_move = min(islas$pob[v] / islas$pob[i], 1)
  if(p_move  > u2){
    return(v) # isla destino
  }
  else{
    return(i) # me quedo en la misma isla
  }
}

pasos <- 100000
camino <- numeric(pasos)
camino[1] <- sample(1:10, 1) # isla inicial
for(j in 2:pasos){
  camino[j] <- caminaIsla(camino[j - 1])
}

caminata <- data.frame(pasos = 1:pasos, isla = camino)

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)))

grid.arrange(plot_caminata, plot_dist, ncol = 1, heights = c(4, 2))
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

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 de cada valor de \(\theta\) simplemente contando el número relativo de veces que el vendedor visitó dicha isla.

t <- c(1:10, 20, 50, 100, 200, 1000, 5000)

plots_list <- lapply(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))
})

args.list <- c(plots_list,list(nrow=4,ncol=4))
do.call(grid.arrange, args.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 de ubica en \(\theta_{actual}\) 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.

  1. 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:

\[p_{mover}=min\bigg( \frac{P(\theta_{propuesta})}{P(\theta_{actual})},1\bigg)\]

Notemos que la distribución objetivo \(P(\theta)\) no necesita estar normalizada, esto es porque lo que nos interesa es el cociente \(P(\theta_{propuesta})/P(\theta_{actual})\).

  1. 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 \(p_{mover}\) entonces hacemos el movimiento.

Entonces, para utilizar el algoritmo necesitamos ser capaces de:

  • Generar un valor de la distribución propuesta (para crear \(\theta_{propuesta}\)).

  • Evaluar la distribución objetivo en cualquier valor propuesto (para calcular \(P(\theta_{propuesta})/P(\theta_{actual})\)).

  • Generar un valor uniforme (para movernos con probabilidad \(p_{mover}\))

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.

library(reshape2)
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 <- ldply(t, function(i){w %*% (T %^% i)})

expT <- expT %>%
  mutate(t = t) %>%
  gather(theta, P, -t) 

ggplot(expT, aes(x = as.numeric(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.

inicioP <- function(i){
  w <- rep(0, 10)
  w[i] <- 1
  t <- c(1, 10, 50, 100)
  expT <- ldply(t, function(i){w %*% (T %^% i)})
  expT <- expT %>%
    mutate(t = t, inicio = i) %>%
    gather(theta, P, -t, -inicio) 
}

expT <- ldply(c(1, 3, 5, 9), inicioP)
ggplot(expT, 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))

Caso general

En la sección anterior implementamos el algoritmo de Metrópolis en un caso sencillo: las posiciones eran discretas, en una dimensión y la propuesta era únicamente mover a la izquierda o a la derecha. El algoritmo general aplica para valores continuos, en cualquier número de dimensiones y con distribuciones propuesta más generales. Lo esencial del método no cambia para el caso general, esto es:

  1. Tenemos una distribución objetivo \(p(\theta)\) de la cual buscamos generar muestras. Debemos ser capaces de calcular el valor de \(p(\theta)\) para cualquier valor candidato \(\theta\). La distribución objetivo \(p(\theta)\) no tiene que estar normalizada, típicamente \(p(\theta)\) es la distribución posterior de \(\theta\) no normalizada, es decir, es el producto de la verosimilitud y la inicial.

  2. La muestra de la distribución objetivo se genera mediante una caminata aleatoria a través del espacio de parámetros. La caminata inicia en un lugar arbitrario (definido por el usuario). El punto inicial debe ser tal que \(p(\theta)>0\). La caminata avanza en cada tiempo proponiendo un movimiento a una nueva posición y después decidiendo si se acepta o no el valor propuesto. Las distribuciones propuesta pueden tener muchas formas, el objetivo es que la distribución propuesta explore el espacio de parámetros de manera eficiente.

  3. Una vez que tenemos un valor propuesto decidimos si aceptar calculando:

\[p_{mover}=min\bigg( \frac{P(\theta_{propuesta})}{P(\theta_{actual})},1\bigg)\]

Y al final obtenemos valores representativos de la distribución objetivo \(\{\theta_1,...,\theta_n\}\)

Es importante recordar que debemos excluir las primeras observaciones pues estas siguen bajo la influencia del valor inicial.

Retomemos el problema de inferencia Bayesiana y veamos como usar el algoritmo de Metrópolis cuando la distribución objetivo es la distribución posterior.

Ejemplo: Bernoulli

Retomemos el ejemplo del experimento Bernoulli, iniciamos con una función de distribución que describa nuestro conocimiento inicial y tal que podamos calcular \(p(\theta)\) con facilidad. En este caso elegimos una densidad beta y podemos usar función dbeta de R:

# p(theta) con theta = 0.4, a = 2, b = 2
dbeta(0.4, 2, 2)
## [1] 1.44
# Definimos la distribución inicial
prior <- function(a = 1, b = 1){
  function(theta) dbeta(theta, a, b)
}

También necesitamos especificar la función de verosimilitud, en nuestro caso tenemos repeticiones de un experimento Bernoulli por lo que: \[\mathcal{L}(\theta) \propto \theta^{11}(1-\theta)^{19}\]

y en R:

# Verosimilitid binomial
likeBern <- function(z, N){
  function(theta){
    theta ^ z * (1 - theta) ^ (N - z)
  }
}

Por tanto la distribución posterior \(p(\theta|x)\) es, por la regla de Bayes, proporcional a \(p(x|\theta)p(\theta)\). Usamos este producto como la distribución objetivo en el algoritmo de Metrópolis.

# posterior no normalizada
postRelProb <- function(theta){
  mi_like(theta) * mi_prior(theta)
}

Implementemos el algoritmo con una inicial Beta(1,1) (uniforme) y observaciones \(z = \sum{x_i}=11\) y \(N = 14\), es decir lanzamos 14 volados de los cuales 11 resultan en águila.

# Datos observados
N = 14
z = 11

# Defino mi inicial y la verosimilitud
mi_prior <- prior() # inicial uniforme
mi_like <- likeBern(z, N) # verosimilitud de los datos observados

# para cada paso decidimos el movimiento de acuerdo a la siguiente función
caminaAleat <- function(theta){ # theta: valor actual
  salto_prop <- rnorm(1, 0, sd = 0.1) # salto propuesto
  theta_prop <- theta + salto_prop # theta propuesta
  if(theta_prop < 0 | theta_prop > 1){ # si el salto implica salir del dominio
    return(theta)
  }
  u <- runif(1) 
  p_move = min(postRelProb(theta_prop) / postRelProb(theta), 1) # prob mover
  if(p_move  > u){
    return(theta_prop) # aceptar valor propuesto
  }
  else{
    return(theta) # rechazar
  }
}

set.seed(47405)

pasos <- 6000
camino <- numeric(pasos) # vector que guardará las simulaciones
camino[1] <- 0.1 # valor inicial

# Generamos la caminata aleatoria
for (j in 2:pasos){
  camino[j] <- caminaAleat(camino[j - 1])
}

caminata <- data.frame(pasos = 1:pasos, theta = camino)

ggplot(caminata[1:3000, ], aes(x = pasos, y = theta)) +
  geom_point(size = 0.8) +
  geom_path(alpha = 0.5) +
  scale_y_continuous(expression(theta), limits = c(0, 1)) +
  scale_x_continuous("Tiempo") +
  geom_vline(xintercept = 600, color = "red", alpha = 0.5)

# excluímos las primeras observaciones (etapa de calentamiento)
caminata_f <- filter(caminata, pasos > 600)
ggplot(caminata_f, aes(x = theta)) +
  geom_density(adjust = 2) +
  labs(title = "Distribución posterior", 
       y = expression(p(theta)), 
       x = expression(theta)) + 
  stat_function(fun = mi_prior, color = "orange") + # inicial
  xlim(0, 1)

Si la distribución objetivo es muy dispersa y la distribución propuesta muy estrecha, entonces se necesitarán muchos pasos para que la caminata aleatoria cubra la distribución con una muestra representativa.

Por otra parte, si la distribución propuesta es muy dispersa podemos caer en rechazar demasiados valores propuestos. Imaginemos que \(\theta_{actual}\) se ubica en una zona de densidad alta, entonces cuando los valores propuestos están lejos del valor actual se ubicarán en zonas de menor densidad y \(p(\theta_{propuesta})/p(\theta_{actual})\) tenderá a ser chico y el movimiento propuesto será aceptado en pocas ocasiones.

De la muestra de valores de \(p(\theta|x)\) obtenidos usando el algoritmo de Metrópolis podemos estimar aspectos de la verdadera distribución \(p(\theta|x)\). Por ejemplo, para resumir la tendencia central es fácil calcular la media y la mediana.

mean(caminata_f$theta)
## [1] 0.748
sd(caminata_f$theta)
## [1] 0.109

En el caso de predicción:

sims_y <- rbinom(nrow(caminata_f), size = 1, prob = caminata_f$theta)
mean(sims_y) # p(y = 1 | x) probabilidad predictiva
## [1] 0.75

Inferencia de dos proporciones binomiales

Consideramos la situación en la que nos interesa estudiar dos proporciones \(\theta_1\) y \(\theta_2\) correspondientes a dos grupos. Queremos determinar que debemos creer de estas proporciones tras observar datos provenientes de ambos grupos.

En el marco Bayesiano comenzamos definiendo nuestras creencias iniciales. En este caso nuestras creencias describen combinaciones de parámetros, es decir debemos especificar \(p(\theta_1, \theta_2)\) para todas las combinaciones \(\theta_1, \theta_2\).

Un caso sencillo es asumir que nuestras creencias de \(\theta_1\) son independientes de nuestras creencias de \(\theta_2\). Esto implica:

\[p(\theta_1, \theta_2) = p(\theta_1)p(\theta_2)\] para todo valor de \(\theta_1\) y de \(\theta_2\) y donde \(p(\theta_1)\) y \(p(\theta_2)\) corresponden a las distribuciones marginales. Las creencias de los dos parámetros no tienen porque ser independientes, por ejemplo, puedo creer que dos monedas acuñadas en la misma casa de moneda tienen un sesgo similar, en este caso las manipulaciones matemáticas son más complicadas pero es posible describir estas creencias mediante una distribución bivariada.

Además de las creencias iniciales tenemos datos observados. Suponemos que los lanzamientos dentro de cada grupo son independientes y que los lanzamientos entre los grupos también lo son. Es importante recalcar que siempre suponemos independencia en los datos, sin importar nuestros supuestos de independencia en las creencias.

Para el primer grupo observamos la secuencia \(\{x_{1,1},...,x_{1,N_1}\}\) que contiene \(z_1\) águilas, y en el otro grupo observamos la sucesión \(\{x_{2,1},...,x_{2,N_2}\}\) que contiene \(z_2\) águilas. Es decir, \[z_1 = \sum_{i=1}^{N_1}x_{1,i}\] Por simplicidad denotamos los datos por \(x =\{x_{1,1},...,x_{1,N_1}x_{2,1},...,x_{2,N_2}\}\)

Debido a la independencia de los lanzamientos tenemos:

\[p(x|\theta_1,\theta_2)=\prod_{i=1}^{N_1}p(x_{1,i}|\theta_1,\theta_2)\cdot \prod_{i=1}^{N_2}p(x_{2,i}|\theta_1,\theta_2)\] \[= \theta_1^{z_1}(1-\theta_1)^{N_1-z_1}\cdot \theta_2^{z_2}(1-\theta_2)^{N_2-z_2}\]

Usamos la regla de Bayes para calcular la distribución posterior:

\[p(\theta_1,\theta_2|x)=p(x|\theta_1,\theta_2)p(\theta_1, \theta_2) / p(x)\] \[= \frac{p(x|\theta_1,\theta_2)p(\theta_1, \theta_2)} { \int\int p(x|\theta_1,\theta_2)p(\theta_1, \theta_2)d\theta_1d\theta_2}\]

Distribuciones conjugadas

Siguiendo el caso de la familia Beta-Bernoulli que estudiamos en el caso de una proporción y suponiendo independencia en las creencias, \(p(\theta_1, \theta_2) = p(\theta_1)p(\theta_2)\). Escribimos la densidad inicial como el producto de dos densidades beta, donde \(\theta_1\) se distribuye \(Beta(a_1, b_1)\) y \(\theta_2\) se distribuye \(Beta(a_2, b_2)\).

\[p(\theta_1, \theta_2)=\frac{\theta_1^{a_1-1}(1-\theta)^{b_1-1}} {B(a_1, b_1)}\frac{\cdot \theta_2^{a_2-1}(1-\theta)^{b_2-1}}{B(a_2, b_2)}\]

donde \(B(a, b) = \Gamma(a)\Gamma(b) / \Gamma(a+b)\).

grid <- expand.grid(x = seq(0.01, 1, 0.01), y = seq(0.01, 1, 0.01))
grid <- grid %>% 
  mutate(z = dbeta(x, 3, 3) * dbeta(y, 3, 3))

ggplot(grid, aes(x = x, y = y, z = z)) + 
  # geom_tile(aes(fill = z)) +
  stat_contour(binwidth = 0.6, aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
  coord_fixed()

La verosimilitud se pude graficar como,

# z_1=5, N_1=7, z_2=2, N_2=7, a_1=a_2=b_1=b_2=3
grid <- grid %>% 
  mutate(z = dbeta(x, 6, 3) * dbeta(y, 3, 6))

ggplot(grid, aes(x = x, y = y, z = z)) + 
  # geom_tile(aes(fill = z)) +
  stat_contour(binwidth = 0.6, aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
  coord_fixed()

la posterior se escribe:

\[p(\theta_1, \theta_2|x)=\frac{\theta_1^{z_1+a_1-1}(1-\theta)^{b_1-1}\cdot \theta_2^{z_2+a_2-1}(1-\theta)^{b_2-1}}{p(x)\cdot B(a_1, b_1) \cdot B(a_2, b_2)}\]

# z_1=5, N_1=7, z_2=2, N_2=7
grid <- grid %>% 
  mutate(z = dbeta(x, 8, 5) * dbeta(y, 5, 8))
  
ggplot(grid, aes(x = x, y = y, z = z)) + 
  # geom_tile(aes(fill = z)) +
  stat_contour(binwidth = 0.6, aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
  coord_fixed()

Resumiendo, cuando la inicial es el producto de distribuciones beta independientes, la posterior también es el producto de distribuciones beta independientes.

Metrópolis

Al igual que en el caso univariado escribimos la verosimilitud y la distribución inicial:

# Verosimilitid binomial
likeBern2 <- function(z_1, N_1, z_2, N_2){
  function(theta){
    theta[1] ^ z_1 * (1 - theta[1]) ^ (N_1 - z_1) *
    theta[2] ^ z_2 * (1 - theta[2]) ^ (N_2 - z_2)
  }
}

prior2 <- function(a_1 = 3, b_1 = 3, a_2 = 3, b_2 = 3){
  function(theta) dbeta(theta[1], a_1 , b_1) * dbeta(theta[2], a_2 , b_2)
}

# posterior no normalizada
postRelProb2 <- function(theta){
  mi_like(theta) * mi_prior(theta)
}

Implementemos el algoritmo con una iniciales Beta(3,3) y observaciones \(z_1=5\), \(N = 7\), \(z_1=2\) y \(N = 7\), es decir lanzamos 14 volados 7 de cada moneda.

z_1=5; N_1=7; z_2=2; N_2=7; a_1=3; a_2=3; b_1=3; b_2=3
# Datos observados
N = 14
z = 11

# Defino mi inicial y la verosimilitud
mi_prior <- prior2() # inicial uniforme
mi_like <- likeBern2(z_1 = z_1, N_1 = N_1, z_2 = z_2, N_2 = N_2) # verosimilitud

Para proponer saltos usaremos una distribución normal bivariada. El movimiento se acepta de manera definitiva si la posterior es más densa que la posición actual y se acepta de manera probabilista si la posición actual es más densa.

# para cada paso decidimos el movimiento de acuerdo a la siguiente función
caminaAleat2 <- function(theta){ # theta: valor actual
    salto_prop <- mvrnorm(n = 1 , mu = rep(0, 2), 
    Sigma = matrix(c(0.1, 0, 0, 0.1), byrow = TRUE, nrow = 2)) # salto propuesto
  theta_prop <- theta + salto_prop # theta propuesta
  if(all(theta_prop < 0) | all(theta_prop > 1)){ # salir del dominio
    return(theta)
  }
  u <- runif(1) 
  p_move = min(postRelProb2(theta_prop) / postRelProb2(theta), 1) # prob mover
  if(p_move  > u){
    return(theta_prop) # aceptar valor propuesto
  }
  else{
    return(theta) # rechazar
  }
}

set.seed(47405)

pasos <- 6000
camino <- matrix(0, nrow = pasos, ncol = 2) # vector que guardará las simulaciones
camino[1, ] <- c(0.1, 0.8) # valor inicial

# Generamos la caminata aleatoria
for (j in 2:pasos){
  camino[j, ] <- caminaAleat2(camino[j - 1, ])
}

caminata <- data.frame(pasos = 1:pasos, theta_1 = camino[, 1], 
                       theta_2 = camino[, 2])

ggplot(caminata, aes(x = theta_1, y = theta_2)) +
  geom_point(size = 0.8) +
  geom_path(alpha = 0.3) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
  coord_fixed()

caminata_m <- caminata %>%
  gather(parametro, val, theta_1, theta_2) %>%
  arrange(pasos)

ggplot(caminata_m[1:2000, ], aes (x = pasos, y = val)) +
  geom_path(alpha = 0.5) +
  facet_wrap(~parametro, ncol = 1) +
  scale_y_continuous("", limits = c(0, 1))

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. Por otra parte, el muestredor de Gibbs no necesita de una distribución propuesta.

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: \[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 Morkov vía Monte Carlo. La diferencia entre Gibbs y Metrópolis radica en como se deciden los pasos.

En el caso de Gibbs, en cada punto de la caminata se selecciona uno de los componentes del vector de parámetros (típicamente se cicla en orden). Supongamos que se selecciona el parámetro \(\theta_k\), entonces obtenemos un nuevo valor para este parámetro generando una simulación de la distribución condicional \[p(\theta_k,|\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p,x)\]

El nuevo valor \(\theta_k\) junto con los valores que aun no cambian \(\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p\) constituyen la nueva posición en la caminata aleatoria.

Seleccionamos una nueva componente 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 si podemos determinar todas las distribuciones condicionales y simular de ellas.

Ejemplificaremos el muestreador de Gibbs con el ejemplo de las proporciones, a pesar de no ser necesario en este caso.

Comenzamos identificando las distribuciones condicionales posteriores para cada parámetro:

\[p(\theta_1|\theta_2,x) = p(\theta_1,\theta_2|x) / p(\theta_2|x)\] \[= \frac{p(\theta_1,\theta_2|x)} {\int p(\theta_1,\theta_2|x) d\theta_1}\]

Usando iniciales \(beta(a_1, b_1)\) y \(beta(a_2,b_2)\), obtenemos:

\[p(\theta_1|\theta_2,x) = \frac{beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}{\int beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2) d\theta_1}\] \[= \frac{beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}{beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}\] \[=beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1)\]

Debido a que la posterior es el producto de dos distribuciones Beta independientes es claro que \(p(\theta_1|\theta_2,x)=p(\theta_1|x)\).

Una vez que determinamos las distribuciones condicionales, simplemente hay que encontrar una manera de obtener muestras de estas, en R podemos usar la función \(rbeta\).

pasos <- 12000
camino <- matrix(0, nrow = pasos, ncol = 2) # vector que guardará las simulaciones
camino[1, 1] <- 0.1 # valor inicial
camino[1, 2] <- 0.1

# Generamos la caminata aleatoria
for (j in 2:pasos){
  if(j %% 2){
    camino[j, 1] <- rbeta(1, z_1 + a_1, N_1 - z_1 + b_1)
    camino[j, 2] <- camino[j - 1, 2]
  }
  else{
    camino[j, 2] <- rbeta(1, z_2 + a_2, N_2 - z_2 + b_2)
    camino[j, 1] <- camino[j - 1, 1]
  }
}

caminata <- data.frame(pasos = 1:pasos, theta_1 = camino[, 1], 
  theta_2 = camino[, 2])

ggplot(caminata[1:2000, ], aes(x = theta_1, y = theta_2)) +
    geom_point(size = 0.8) +
    geom_path(alpha = 0.5) +
    scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
    scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
    coord_fixed()

caminata_g <- filter(caminata, pasos %% 2 == 0) %>%
  gather(parametro, val, theta_1, theta_2) %>%
  mutate(pasos = rep(1:6000, 2)) %>%
  arrange(pasos)

ggplot(caminata_g[1:2000, ], aes(x = pasos, y = val)) +
  geom_path(alpha = 0.3) +
  facet_wrap(~parametro, ncol = 1) +
  scale_y_continuous("", limits = c(0, 1))

Si comparamos los resultados del muestreador de Gibbs con los de Metrópolis notamos que las estimaciones son muy cercanas

# Metropolis
caminata_m %>% 
  filter(pasos > 1000) %>% # eliminamos el calentamiento
  group_by(parametro) %>%
  summarise(
    media = mean(val),
    mediana = median(val),
    std = sd(val)
    )
## Source: local data frame [2 x 4]
## 
##   parametro media mediana   std
## 1   theta_1 0.605   0.606 0.124
## 2   theta_2 0.392   0.388 0.136
# Gibbs
caminata_g %>% 
  filter(pasos > 1000) %>%
  group_by(parametro) %>%
  summarise(
    media = mean(val),
    mediana = median(val),
    std = sd(val)
    )
## Source: local data frame [2 x 4]
## 
##   parametro media mediana   std
## 1   theta_1 0.619   0.630 0.129
## 2   theta_2 0.384   0.378 0.131

También podemos comparar los sesgos de las dos monedas, esta es una pregunta más interesante.

caminata <- caminata %>%
  mutate(dif = theta_1 - theta_2)

ggplot(caminata, aes(x = dif)) + 
  geom_histogram(fill = "gray") + 
  geom_vline(xintercept = 0, alpha = 0.8, color = "red")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

La principal ventaja del muestreador de Gibbs sobre el algoritmo de Metrópolis es que no hay necesidad de seleccionar una distribución propuesta y no hay que lidiar con lo ineficiente de rechazar valores. A cambio, debemos ser capaces de derivar las probabilidades condicionales de cada parámetro y de generar muestras de estas.

Ejemplo: Normal

Retomemos el caso de observaciones normales, supongamos que tengo una muestra \(x_1,...,x_N\) de observacionesindependientes e identicamente distribuidas, con \(x_i \sim N(\mu, \sigma^2)\), veremos el caso de media desconocida, varianza desconocida y de ambas desconocidas.

Normal con media desconocida. Supongamos que \(\sigma^2\) es conocida, por lo que nuestro parámetro de interés es únicamente \(\mu\) entonces si describo mi conocimiento inicial de \(\mu\) a través de una distribución normal: \[\mu \sim N(m, \tau^2)\] resulta en una distribución posterior: \[\mu|x \sim N\bigg(\frac{\sigma^2}{\sigma^2 + N\tau^2}m + \frac{N\tau^2}{\sigma^2 + N \tau^2}\bar{x}, \frac{\sigma^2 \tau^2}{\sigma^2 + N\tau^2}\bigg)\]

Normal con varianza desconocida. Supongamos que \(\mu\) es conocida, por lo que nuestro parámetro de interés es únicamente \(\sigma^2\). En este caso una distribución conveniente para describir nuestro conocimiento inicial es la distribución Gamma Inversa.

La distribución Gamma Inversa es una distribución continua con dos parámetros y que toma valores en los positivos. Como su nombre lo indica, esta distribución corresponde al recírpoco de una variable cuya distribución es Gamma, recordemos que si \(x\sim Gamma(\alpha, \beta)\) entonces:

\[p(x)=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}x^{\alpha-1}e^{-x/\beta}\]

donde \(x>0\). Ahora si \(y\) es la variable aleatoria recírpoco de \(x\) entonces:

\[p(y)=\frac{\beta^\alpha}{\Gamma(\alpha)}y^{-\alpha - 1} exp{-\beta/y}\]

con media \[\frac{\beta}{\alpha-1}\] y varianza \[\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}.\]

Debido a la relación entre las distribuciones Gamma y Gamma Inversa, podemos utilizar la función rgamma de R para generar valores con distribución gamma inversa.

# 1. simulamos valores porvenientes de una distribución gamma
x_gamma <- rgamma(2000, shape = 5, rate = 1)
# 2. invertimos los valores simulados
x_igamma <- 1 / x_gamma

# También podemos usar las funciones de MCMCpack
library(MCMCpack)
## Loading required package: coda
## Loading required package: lattice
## 
## Attaching package: 'coda'
## 
## The following object is masked from 'package:arm':
## 
##     traceplot
## 
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2015 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
x_igamma <- data.frame(x_igamma)

ggplot(x_igamma, aes(x = x_igamma)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, fill = "gray") + 
  stat_function(fun = dinvgamma, args = list(shape = 5, scale = 1), 
    color = "red")

Volviendo al problema de inferir acerca del parámetros \(\sigma^2\), si resumimos nuestro conocimiento inicial a través de una distribución Gamma Inversa tenemos \[p(\sigma^2)=\frac{\beta^\alpha}{\Gamma(\alpha)}\frac{1}{(\sigma^2)^{\alpha + 1}} e^{-\beta/\sigma^2}\]

la verosimiltud: \[p(x|\mu, \sigma^2)=\frac{1}{(2\pi\sigma^2)^{N/2}}exp\left(-\frac{1}{2\sigma^2}\sum_{j=1}^{N}(x_j-\mu)^2\right)\]

y calculamos la posterior:

\[p(\sigma^2) \propto p(x|\mu,\sigma^2)p(\sigma^2)\]

obtenemos que \(\sigma^2|x \sim GI(N/2+\alpha, \beta + 1/2 \sum(x_i - \mu)^2)\).

Por tantyo tenemos que la inicial Gamma con verosimilitud Normal es una familia conjugada.

Normal con media y varianza desconocidas. Sea \(\theta=(\mu, \sigma^2)\) especificamos la siguiente inicial para \(\theta\): \[p(\theta) = N(\mu|m, \tau^2)\cdot IG(\sigma^2|\alpha, \beta)\] suponemos hiperparámetros \(m,\tau^2, \alpha, \beta\) conocidos. Entonces, la distribución posterior es: \[ p(\theta|x) \propto p(x|\theta) p(\theta)\] \[= \frac{1}{(\sigma^2)^{N/2}} exp\bigg(-\frac{1}{2\sigma^2}\sum_{i=1}^N (x_i-\mu)^2 \bigg) exp\bigg(-\frac{1}{2\tau^2}(\mu-m)^2)\bigg) \frac{1}{(\sigma^2)^{\alpha +1}} exp\bigg(-\frac{\beta}{\sigma^2}\bigg) \]

en esta última distribución no reconocemos el núcleo de niguna distribución conocida pero si nos concenteramos únicamente en los términos que involucran a \(\mu\) tenemos: \[exp\left(-\frac{1}{2}\left( \mu^2 \left( \frac{N}{\sigma^2} + \frac{1}{\tau^2} \right) - 2\mu\left(\frac{\sum_{i= 1}^n x_i}{\sigma^2} + \frac{m}{\tau^2}\right) \right)\right) \] esta expresión depende de \(\mu\) y \(\sigma^2\), sin embargo condicional a \(\sigma^2\) observamos el núcleo de una distribución normal, \[\mu|\sigma^2,x \sim N\left(\frac{n\tau^2}{n\tau^2 + \sigma^2}\bar{x} + \frac{\sigma^2}{N\tau^2 + \sigma^2}m, \frac{\tau^2\sigma^2}{n\tau^2 + \sigma^2} \right)\] Si nos fijamos únicamente en los tárminos que involucran a \(\sigma^2\) tenemos: \[\frac{1}{(\sigma^2)^{N/2+\alpha+1}}exp\left(- \frac{1}{\sigma^2} \left(\sum_{i=1}^N \frac{(x_i-\mu)^2}{2} + \beta \right) \right)\] y tenemos \[\sigma^2|\mu,x \sim GI\left(\frac{N}{2} + \alpha, \sum_{i=1}^n \frac{(x_i-\mu)^2}{2} + \beta \right)\] Obtenemos así las densidades condicionales completas \(p(\mu|\sigma^2, x)\) y \(p(\sigma^2|\mu, x)\) cuyas distribuciones conocemos y de las cuales podemos simular.

Implementaremos un muestreador de Gibbs.

Comenzamos definiendo las distrbuciones iniciales:

  • \(\mu \sim N(1.5, 16)\), esto es \(m = 1.5\) y \(\tau^2 = 16\).

  • \(\sigma^2 \sim GI(3, 3)\), esto es \(\alpha = \beta = 3\).

Ahora supongamos que observamos 20 realizaciones provenientes de la distribución de interés:

N <- 50 # Observamos 20 realizaciones
set.seed(122)
x <- rnorm(N, 2, 2) 
x
##  [1]  4.6214  0.2483  2.3990  2.9319 -1.6041  4.8975  2.5977  2.7236
##  [9] -0.0139  1.4860  1.7357  0.3167  2.5485 -2.9252 -2.3068  4.3184
## [17]  3.3795  3.7605  0.1133  3.4381  0.9243  0.9547 -0.1058  2.2030
## [25]  5.7270  1.9608 -0.1566  2.3452  3.0661  5.9045  4.8227  3.2027
## [33]  0.1720  5.1609  1.0602  5.2037  2.7455  2.0678  2.2082 -2.0367
## [41]  0.6840 -1.2170 -0.6882  1.6067  4.7513  2.3466  1.1214 -1.1906
## [49]  0.5114  5.6304

Escribimos el muestreador de Gibbs.

m <- 1.5; tau2 <- 16; alpha <- 3; beta <- 3 # parámetros de iniciales

pasos <- 20000
camino <- matrix(0, nrow = pasos + 1, ncol = 2) # vector guardará las simulaciones
camino[1, 1] <- 0 # valor inicial media

# Generamos la caminata aleatoria
for (j in 2:(pasos + 1)){
  # sigma^2
  mu <- camino[j - 1, 1]
  a <- N / 2 + alpha
  b <- sum((x  - mu) ^ 2) / 2 + beta
  camino[j - 1, 2] <- 1/rgamma(1, shape = a, rate = b) # Actualizar sigma2
  
  # mu
  sigma2 <- camino[j - 1, 2]
  media <- (N * tau2 * mean(x) + sigma2 * m) / (N * tau2 + sigma2)
  var <- sigma2 * tau2 / (N * tau2 + sigma2)
  camino[j, 1] <- rnorm(1, media, sd = sqrt(var)) # actualizar mu
}

caminata <- data.frame(pasos = 1:pasos, mu = camino[1:pasos, 1], 
  sigma2 = camino[1:pasos, 2])

caminata_g <- caminata %>%
  gather(parametro, val, mu, sigma2) %>%
  arrange(pasos)

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("")

ggplot(filter(caminata_g, pasos > 5000), aes(x = val)) +
  geom_histogram(fill = "gray") +
  facet_wrap(~parametro, ncol = 1, scales = "free") 
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

Algunos resúmenes de la posterior:

caminata_g %>%
  filter(pasos > 1000) %>% # eliminamos la etapa de calentamiento
  group_by(parametro) %>%
  summarise(
    mean(val), 
    sd(val), 
    median(val)
    )
## Source: local data frame [2 x 4]
## 
##   parametro mean(val) sd(val) median(val)
## 1        mu      1.91   0.305        1.91
## 2    sigma2      4.74   0.933        4.62

Predicción. Para predecir el valor de una realización futura \(y\) recordemos que:

\[p(y) =\int p(y|\theta)p(\theta|x)d\theta\]

Por tanto podemos aproximar la distribución predictiva posterior como:

caminata_f <- filter(caminata, pasos > 5000)

caminata_f$y_sims <- rnorm(1:nrow(caminata_f), caminata_f$mu, caminata_f$sigma2)

ggplot(caminata_f, aes(x = y_sims)) + 
  geom_histogram(fill = "gray") +
  geom_vline(aes(xintercept = mean(y_sims)), color = "red")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

¿Cuál es la probabilidad de que una observación futura sea mayor a 8?

En estadística bayesiana es común parametrizar la distribución Normal en términos de precisión (el inverso de la varianza). Si parametrizamos de esta manera \(\nu = 1/\sigma^2\) podemos repetir el proceso anterior con la diferencia de utilizar la distribución Gamma en lugar de Gamma inversa.

JAGS

  • Instalar JAGS.

  • Instalar los paquetes R2jags y rjags de R.

JAGS (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.

Repitamos nuestro ejemplo usando JAGS, es importante tomar en cuenta que en JAGS la distribución Normal esta parametrizada con la presición.

La distribución inicial y la verosimilitud se especifican en el modelo:

modelo_normal.txt <-
'
model{
  for(i in 1:N){
    x[i] ~ dnorm(mu, nu)
  }
  # iniciales
  nu ~ dgamma(3, 3)
  sigma2 <- 1 / nu
  mu ~ dnorm(1.5, 1 / 16) 
}
'

el ciclo for indica que cada dato observado \(x_i\) proviene de una distribución Normal con media \(\mu\) y varinza \(1 / \nu\) (precisión \(\nu\)). Afuera del ciclo escribimos las distribuciones iniciales, \(\nu \sim Gamma(3, 3)\), esto es \(\sigma^2 \sim GI(3, 3)\) y \(\mu\) se distribuye Normal con media \(\mu = 1.5\) y varianza \(\tau^2 = 16\).

El modelo ya esta especificado, aun debemos indicar los valores de las variables en el modelo, para esto definimos los valores en R y después los mandamos a JAGS.

library(R2jags)
## Loading required package: rjags
## Linked to JAGS 3.4.0
## Loaded modules: basemod,bugs
## 
## Attaching package: 'R2jags'
## 
## The following object is masked from 'package:coda':
## 
##     traceplot
## 
## The following object is masked from 'package:arm':
## 
##     traceplot
cat(modelo_normal.txt, file = 'modelo_normal.bugs')

# valores iniciales para los parámetros, si no se especifican la función jags
# generará valores iniciales
jags.inits <- function(){
  list("mu" = rnorm(1, 0, 20),
    "nu" = 1/runif(0, 100))
}

jags_fit <- jags(
  model.file = "modelo_normal.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, N = N),    # lista con los datos
  parameters.to.save = c("mu", "nu", "sigma2"),  # parámetros por guardar
  n.chains = 1,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000,   # calentamiento de la cadena
  n.thin = 1
  )
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 59
## 
## Initializing model
jags_fit
## Inference for Bugs model at "modelo_normal.bugs", fit using jags,
##  1 chains, each with 10000 iterations (first 1000 discarded)
##  n.sims = 9000 iterations saved
##          mu.vect sd.vect    2.5%     25%     50%     75%  97.5%
## mu         1.914   0.307   1.314   1.709   1.912   2.114   2.52
## nu         0.219   0.042   0.145   0.189   0.216   0.245   0.31
## sigma2     4.734   0.938   3.230   4.074   4.625   5.282   6.88
## deviance 223.464   2.060 221.455 221.997 222.811 224.303 228.78
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2.1 and DIC = 225.6
## DIC is an estimate of expected predictive error (lower deviance is better).
# podemos ver las cadenas
traceplot(jags_fit, varname = c("mu", "sigma2"))

names(jags_fit$BUGSoutput)
##  [1] "n.chains"        "n.iter"          "n.burnin"       
##  [4] "n.thin"          "n.keep"          "n.sims"         
##  [7] "sims.array"      "sims.list"       "sims.matrix"    
## [10] "summary"         "mean"            "sd"             
## [13] "median"          "root.short"      "long.short"     
## [16] "dimension.short" "indexes.short"   "last.values"    
## [19] "program"         "model.file"      "isDIC"          
## [22] "DICbyR"          "pD"              "DIC"
  • Realiza un histograma de la distribución predictiva. Construye un intervalo de 95% de probabilidad para la predicción. Pista: utiliza jags_fit\(BUGSoutput\)sims.matrix.
jags_fit <- jags(
  model.file = "modelo_normal.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = c(NA, x), N = N + 1),    # lista con los datos
  parameters.to.save = c("mu", "nu", "sigma2", "x"),  # parámetros por guardar
  n.chains = 1,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000,   # calentamiento de la cadena
  n.thin = 1
  )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 60
## 
## Initializing model
head(jags_fit$BUGSoutput$summary)
##             mean     sd    2.5%     25%     50%     75%   97.5%
## deviance 223.447 2.0372 221.462 221.981 222.824 224.231 228.854
## mu         1.912 0.3075   1.298   1.712   1.910   2.115   2.520
## nu         0.219 0.0418   0.144   0.189   0.217   0.245   0.309
## sigma2     4.746 0.9447   3.235   4.078   4.617   5.282   6.927
## x[1]       1.930 2.2047  -2.378   0.464   1.914   3.402   6.361
## x[2]       4.621 0.0000   4.621   4.621   4.621   4.621   4.621
mus <- jags_fit$BUGSoutput$sims.matrix[, 2]
sigmas <- sqrt(jags_fit$BUGSoutput$sims.matrix[, 4])
y <- rnorm(length(mus), mus[1], sigmas[1])

mean(y)
## [1] 1.48
sd(y)
## [1] 2.12

Diagnósticos

Para determinar la convergencia de la cadena es conveniente realizar más de dos cadenas, la idea es ver si realmente se ha olvidado el estado inicial viendo si las simulaciones de distintas cadenas convergen a una distribución comín.

jags_fit_pred <- jags(
  model.file = "modelo_normal.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, N = N),    # lista con los datos
  parameters.to.save = c("mu", "nu", "sigma2"),  # parámetros por guardar
  n.chains = 3,   # número de cadenas
  n.iter = 50,    # número de pasos
  n.burnin = 0,   # calentamiento de la cadena
  n.thin = 1
  )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 59
## 
## Initializing model
# podemos ver las cadenas
traceplot(jags_fit, varname = c("mu", "sigma2"))

además de realizar gráficas podemos usar la medida de convergencia \(\hat{R}\) que la función_jags_ calcula por omisión:

jags_fit_pred$BUGSoutput$summary
##             mean     sd   2.5%     25%     50%     75%   97.5%  Rhat n.eff
## deviance 223.461 1.8667 221.45 222.069 222.944 224.376 227.380 1.017   150
## mu         1.921 0.2984   1.34   1.729   1.915   2.131   2.485 0.993   150
## nu         0.218 0.0429   0.14   0.185   0.215   0.251   0.305 1.006   150
## sigma2     4.772 0.9824   3.27   3.979   4.655   5.412   7.127 1.006   150

La medida \(\hat{R}\) se conoce como el factor de reducción potencial de escala este es 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.

Adicionalmente las salidas de la función jags incluyen \(n.eff\) que es el número efectivo de simulaciones, si las simulaciones fueran independientes, \(n.eff\) sería el número total de simulaciones; sin embargo, las simulaciones de MCMC suelen estar correlacionadas y por tanto el número efetivo de simulaciones es menor. Usualmente nos gustaría obtener un número efectivo de al menos \(100\).

Verosimilitud Bernoulli con inicial jerárquica

Ahora estudiaremos el caso en el que los parámetros son dependientes. Podríamos reflejar el caso en que dos monedas provienen de la misma casa de moneda. En este caso tenemos:

  1. Conocimientos iniciales de los posibles valores de los parámetros (sesgos de las monedas).

  2. Tenemos conocimiento inicial de la dependencia de los parámetros por provenir de la misma fábrica.

Cuando observamos lanzamientos de las monedas actualizamos nuestras creencias relativas a los sesgos de las monedas y también actualizamos nuestras creencias acerca de la dependencia de los sesgos.

Suponemos una distribución inicial beta. Recordemos que la distribución beta tienen dos parámetors \(a\) y \(b\):

\[p(\theta)=\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\]

Con el fin de hacer los parámetros más intuitivos los podemos expresar en términos de la media \(\mu\) y el tamaño de muestra \(K\), donde \(\mu\) es la media de nuestro conocimiento inicial y la confianza está reflejada en el tamaño de muestra \(K\). Entonces los parámetros correspondientes en la distribución beta son:

\[a=\mu K, b = (1-\mu)K\]

Ahora introducimos el modelo jerárquico. En lugar de especificar un valor particular para \(\mu\), consideramos que \(\mu\) puede tomar distintos valores (entre 0 y 1), y definimos una distribución de probabilidad sobre esos valores. Podemos pensar que esta distribución describe nuestra incertidumbre acerca de la construcción de la máquina que manufactura las monedas. En nuestro ejemplo, supongamos que la distribución de \(\mu\) es nuevamente beta,

\[p(\mu)=beta(\mu|A_{\mu}, B_{\mu})\]

donde \(A_{\mu}\) y \(B_{\mu}\) se conocen como hiperparámetros y son constantes. En este caso, consideramos que \(\mu\) se ubica típicamente cerca de \(A_{\mu}/(A_{\mu} + B_{\mu})\) y \(K\) se considera constante.

Veamos como aplica la regla de Bayes, si consideramos este como un problema donde tenemos dos parámetros desconocidos \(\mu\) y \(\theta\) tenemos

\[p(\theta, \mu)=\frac{p(x|\theta,\mu)}{p(\theta,\mu)}{p(x)}\]

Hay dos aspectos a considerar en el problema:

  1. La verosimilitud no depende de \(\mu\) por lo que \[p(x|\theta, \mu)=p(x|\theta)\]

  2. La distribución inicial en el espacio de parámetros bivariado se puede factorizar:

\[p(\theta,\mu)=p(\theta|\mu)p(\mu)\]

Por lo tanto \[p(\theta,\mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\] \[=\frac{p(x|\theta,\mu)p(\theta|\mu)p(\mu)}{p(x)}\]

Aproximación por cuadrícula

Si los parámetros e hiperparámetros toman un número finito de valores y no hay muchos de ellos, podemos aproximar la posterior usando aproximación por cuadrícula.

A continuación graficamos las distribuciones correspondientes al caso en que la distribución del hiperparámetro \(\mu\) tiene la forma de una distribución Beta(2, 2), es decir creemos que la media de la acuñadora \(\mu\) es 0.5, pero existe bastante incertidumbre acerca del valor.

La distribución de \(\theta\), esto es la distribución inicial que refleja la dependencia entre \(\theta\) y \(\mu\) se expresa por medio de otra distribución beta:

\[\theta|\mu \approx beta(\mu 100, (1-\mu)100)\]

Esta inicial expresa un alto grado de certeza que una acuñadora con hiperparámetro \(\mu\) genera monedas con sesgo cercano a \(\mu\)

Verosimilitud.

Posterior.

Multiples monedas de una misma fábrica

La sección anterior considera el escenario en que lanzamos una moneda y hacemos inferencia del parámetro de sesgo \(\theta\) y del hiperparámetro \(\mu\). Ahora consideramos recolectar información de múltiples monedas, si cada moneda tiene su propio sesgo \(\theta_j\) entonces debemos estimar un parámetro distinto para cada moneda.

Suponemos que todas las monedas provienen de la misma fábrica, esto implica que tenemos la misma información inicial \(\mu\) para todas las monedas. Suponemos también que cada moneda se acuña de manera independiente, esto es, que condicional al parámetro \(\mu\) los parámetros \(\theta_j\) son independientes en nuestros conocimientos iniciales.

Posterior vía aproximación por cuadrícula

Supongamosque tenemos dos monedas de la misma fábrica. El objetivo es estimar los sesgos \(\theta_1\), \(\theta_2\) de las dos monedas y estimar simultáneamente el parámetro \(\mu\) correspondiente a la casa de moneda que las fabricó.

Inicial.

Verosimilitud.

Posterior.

Muestreador de Gibbs

La sección anterior utilizó un modelo simplificado con el objetivo de poder visualizar el espacio de parámetros. Ahora incluiremos más parámetros con el objetivo de hacer el problema más realista. En los ejemplos anteriores fijamos el grado de dependencia entre \(\mu\) y \(\theta\) de manera arbitraria, a través de \(K\), de tal manera que si \(K\) era grande, los valores \(\theta_j\) individuales se situaban cerca de \(\mu\) y cuando \(K\) era pequeña se permitía más variación.

En situaciones reales no conocemos el valor de \(K\) por adelantado, por lo que dejamos que los datos nos informen acerca de valores creíbles para \(K\). Intuitivamente, cuando la proporción de águilas observadas es similar a lo largo de las monedas, tenemos evidencia de que \(K\) es grande, mientras que si estas proporciones difieren mucho, tenemos evidencia de que \(K\) es pequeña. Debido a que \(K\) pasará de ser una constante a ser un parámetro lo llamaremos \(\kappa\).

La distribución inicial para \(\kappa\) será una Gamma.

modelo_jer.txt <-
'
model{
  for(t in 1:N){
    x[t] ~ dbern(theta[coin[t]])
  }
  for(j in 1:nCoins){
    theta[j] ~ dbeta(a, b)
  }
  a <- mu * kappa
  b <- (1 - mu) * kappa
  # A_mu = 2, B_mu = 2
  mu ~ dbeta(2, 2)
  kappa ~ dgamma(1, 0.1)
}
'

\(\kappa ~Gamma(1, 0.1)\), esta tiene media 10 y desviación estándar 10.

Los datos consisten en 5 monedas, cada una se lanza 5 veces, resultando 4 de ellas en 1 águila y 4 soles y otra en 3 águilas y 2 soles.

cat(modelo_jer.txt, file = 'modelo_jer.bugs')

x <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1)
coin <- c(rep(1, 5), rep(2, 5), rep(3, 5), rep(4, 5), rep(5, 5))

jags.inits <- function(){
  list("mu" = runif(1, 0.1, 0.9),
    "kappa" = runif(1, 5, 20))
}

jags_fit <- jags(
  model.file = "modelo_jer.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, coin = coin, nCoins = 5,  N = 25),    # lista con los datos
  parameters.to.save = c("mu", "kappa", "theta"),  # parámetros por guardar
  n.chains = 3,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000   # calentamiento de la cadena
  )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 65
## 
## Initializing model
traceplot(jags_fit, varname = c("kappa", "mu", "theta"))

jags_fit
## Inference for Bugs model at "modelo_jer.bugs", fit using jags,
##  3 chains, each with 10000 iterations (first 1000 discarded), n.thin = 9
##  n.sims = 3000 iterations saved
##          mu.vect sd.vect   2.5%    25%    50%    75%  97.5% Rhat n.eff
## kappa     14.985  11.497  2.385  6.929 12.025 19.664 44.887    1  3000
## mu         0.330   0.094  0.162  0.263  0.326  0.389  0.525    1  3000
## theta[1]   0.284   0.126  0.062  0.191  0.276  0.365  0.553    1  3000
## theta[2]   0.286   0.126  0.071  0.194  0.277  0.365  0.564    1  3000
## theta[3]   0.285   0.127  0.064  0.193  0.274  0.366  0.557    1  2800
## theta[4]   0.284   0.128  0.069  0.189  0.275  0.367  0.562    1  3000
## theta[5]   0.412   0.141  0.169  0.312  0.399  0.504  0.707    1  1500
## deviance  30.280   2.059 27.450 28.807 29.924 31.319 35.388    1  2800
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 2.1 and DIC = 32.4
## DIC is an estimate of expected predictive error (lower deviance is better).
sims_df <- data.frame(n_sim = 1:jags_fit$BUGSoutput$n.sims,
                      jags_fit$BUGSoutput$sims.matrix) %>% 
  dplyr::select(-deviance) %>%
  gather(parametro, value, -n_sim)

ggplot(filter(sims_df, parametro != "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8) + facet_wrap(~ parametro)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

ggplot(filter(sims_df, parametro == "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8)
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.