¿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 l