In [1]:
%matplotlib inline
import numpy as np
import pandas as pd
import arviz as az
import scipy.stats as stats
from matplotlib import pyplot as plt
from ipywidgets import interact
import ipywidgets as ipyw
az.style.use('arviz-darkgrid')
La estadística trata sobre la recolección, organización, análisis e interpretación de datos, es por ello que la estadística es esencial para el correcto análisis de datos. Existen dos grandes conjuntos de herramientas para analizar datos:
Análisis Exploratorio de Datos (EDA): Consiste en resúmenes numéricos como la media, moda, desviación estándar, rangos intercuartiles, etc (esto se conoce también como estadística descriptiva). Además hace énfasis en el uso de métodos visuales para inspeccionar los datos, como por ejemplo histogramas y gráficos de dispersión.
Estadística Inferencial: Consiste en usar datos para generar enunciados que exceden los propios datos. A veces esto implica realizar predicciones, a veces entender los detalles de algún fenómeno en particular o elegir entre varias explicaciones plausibles.
El foco de este curso está en aprender a realizar inferencias Bayesiana y luego usar métodos del análisis exploratorio de datos para resumir, interpretar, evaluar y comunicar los resultados de tales inferencias.
Muchos de los cursos y libros sobre estadística, al menos los dirigidos a no-estadísticos, enseñan una serie de recetas que más o menos tienen la siguiente forma.
La principal meta de estos cursos es aprender a usar la lata adecuada y algo sobre el emplatado y la presentación del menú. Personalmente nunca me gustó esta aproximación, la razón principal es que el resultado más común de esta forma de enseñanza es un montón de gente confundida incapaces de percibir o entender conceptualmente la unidad de los diferentes métodos enseñados. Por lo tanto nosotros tomaremos una aproximación diferente. También aprenderemos recetas, pero intentaremos que los platos tengan un sabor más casero y menos enlatado, aprenderemos a mezclar ingredientes frescos que se acomoden a diferentes situaciones gastronómicas.
Este enfoque es posible por dos razones:
Probability theory is nothing but common sense reduced to calculation. -Pierre-Simon Laplace
En este curso aprenderemos sobre una forma de hacer estadística llamada usualmente estadística Bayesiana. El nombre se debe a Thomas Bayes (1702-1761) un ministro presbiteriano, y matemático aficionado, quien derivó por primera vez lo que ahora conocemos como el teorema de Bayes, el cual fue publicado (postumanente) en 1763. Sin embargo una de las primeras personas en realmente desarrollar métodos Bayesianos, fue Pierre-Simon Laplace (1749-1827), por lo que tal vez sería un poco más correcto hablar de Estadística Laplaciana y no Bayesiana.
Existe otro paradigma estadístico llamado estadística clásica o frecuentista. Si ustedes han tenido un curso de estadística (ya sea en el grado o posgrado) es casi seguro que dicho curso fue sobre métodos frecuentistas (aun cuando esto no haya sido explicitado). Es interesante notar que mientras los orígenes de las estadística Bayesiana se remontan al siglo XVIII. Los métodos "clásicos" (o frecuentistas) fueron desarrollados principalmente durante el siglo XX! De hecho una de las motivaciones para desarrollar métodos frecuentistas fue un sentimiento e ideología anti-bayesiano. A lo largo del curso nos centraremos en los métodos Bayesianos.
Hay dos ideas centrales que hacen que un método sea Bayesiano:
En el universo Bayesiano las cantidades conocidas son consideradas fijas y usualmente les llamamos datos. Por el contrario toda cantidad desconocida es considerada como una variable aleatoria y modelada usando una distribución de probabilidad.
Las probabiliades son números entre 0 y 1 (incluyendo ambos extremos). En estadística Bayesiana las probabilidades son usadas para cuantificar la confianza que tenemos en que un evento ocurra. Desde este punto de vista es totalmente razonable preguntar cual es la probabilidad de que la masa de Saturno sea $5 \times 10^{26}$ kg, o hablar sobre la probabilidad de lluvia durante el 25 de Mayo de 1810, o la probabilidad de que mañana amanezca.
La lógica aristotélica permite razonar de forma correcta cuando los enunciados son verdaderos o falsos (cuando hay certezas). A fin de razonar, de forma correcta, en presencia de incertidumbre es necesario extender la lógica aristotélica. En 1946, Richard Cox demostró no solo que tal extensión es posible si no que esa extensión se corresponde con la ya conocida teoría de probabilidad.
Si no estamos seguros de la factibilidad de un evento, entonces es matemáticamente razonable asignar un valor entre 0 y 1, de acuerdo al grado de confianza que tenemos de que ocurra dicho evento. Los extremos son los casos especiales de absoluta certeza $p(A=0)= falso$ y $p(A=1) = verdadero$.
Según esta interpretación, llamada Bayesiana, las probabilidades son una medida de la incertidumbre con la que conocemos una cantidad y no necesariamente una propiedad de la naturaleza, distintas personas podrán asignar distintas probabilidades a un mismo evento. Por ello se suele decir que la estadística Bayesiana es subjetiva (¡y no como un cumplido!). Sin embargo, este uso de las probabilidades simplemente refleja un hecho trivial, distintas personas poseen diferente información (la cual bien podría ser rotulada de objetiva) por lo que no necesariamente estarán de acuerdo en el grado de factibilidad de un evento. Una persona que sabe que una moneda está sesgada hacia cara, maneja información distinta de otra que asume que una moneda tiene igual posibilidad de caer cara o seca. De todas formas, si ambas personas hacen el experimento de arrojar esa moneda al aire varias veces, la probabilidad que cada uno asigna al evento cara irá convergiendo a un mismo valor, al menos si se comportan de forma racional. Más adelante realizaremos tal experimento de forma computacional.
Como se verá en el resto del curso la estadística Bayesiana (al menos en su sentido moderno dentro de las ciencias exactas, naturales y sociales) es tan subjetiva u objetiva como los métodos estadísticos no-Bayesianos u otras formas de modelado usadas en ciencia.
Para empezar a ponernos un poco más concretos veamos el teorema de Bayes que, junto con la noción de probabilidad que acabamos de discutir, es la base de la estadística Bayesiana.
Si una variable puede tomar más de un valor $\{x_1, x_2, \cdots, x_n\}$ es posible usar una distribución de probabilidad $p(x)$ para modelar la plausibilidad de cada uno de esos valores. Al trabajar con probabilidades es común definir los siguientes conceptos (resumidos en la siguiente figura).
En estadística Bayesiana la cantidad central a determinar es:
\begin{equation} p(\theta \mid y) \end{equation}Donde $\theta$ es un conjunto de parámetros que se relacionan con nuestro problema e $y$ representa las observaciones. En algún momento de la historia la estadística Bayesiana recibió el nombre de probabilidad inversa, ya que infiere a partir de las observaciones los parámetros. Es importante notar que $\theta$ representa una distribución de valores plausibles y no a una estimación puntual, es decir no es un único número.
El teorema de Bayes es una consecuencia directa de la regla del producto, veamos.
\begin{equation} p(\theta, y) = p(\theta \mid y) p(y) \\ p(\theta, y) = p(y \mid \theta) p(\theta) \end{equation}Dado que los dos términos a la derecha de la igualdad son iguales entre si podemos escribir que:
\begin{equation} p(\theta \mid y) p(y) = p(y \mid \theta) p(\theta) \end{equation}Reordenando llegamos al Teorema de Bayes!
\begin{equation} p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y)} \end{equation}El cual también suele ser escrito de la siguiente forma:
\begin{equation} p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{\int_{\theta} p(y \mid \theta) p(\theta) \text{d}\theta} \end{equation}El denominador es una integral sobre todos los valores de $\theta$ definidos por el a priori. La integral es reemplazada por una sumatoria en el caso que estemos hablando de valores discretos y no continuos.
Cada término del teorema de Bayes tiene un nombre específico:
El a priori es la forma de introducir conocimiento previo sobre los valores que pueden tomar los parámetros. A veces cuando no sabemos demasiado se suelen usar a prioris que asignan igual probabilidad a todos los valores de los parámetros, otras veces se puede elegir a prioris que restrijan los valores de los parámetros a rangos razonables, algo que se conoce como regularización, por ejemplo solo valores positivos. Muchas veces contamos con información mucho más precisa como medidas experimentales previas o límites impuesto por alguna teoría.
El likelihood es la forma de incluir nuestros datos en el análisis. Es una expresión matemática que especifica la plausibilidad de los datos. El likelihood es central tanto en estadística Bayesiana como en estadística no-Bayesiana. A medida que la cantidad de datos aumenta el likelihood tiene cada vez más peso en los resultados, esto explica el porqué a veces los resultados de la estadística Bayesiana y frecuentista coinciden cuando la muestra es grande.
El a posteriori es la distribución de probabilidad para los parámetros. Es la consecuencia lógica de haber usado un conjunto de datos, un likelihood y un a priori. Se lo suele pensar como la versión actualizada del a priori. De hecho un a posteriori puede ser un a priori de un análisis a futuro.
La likelihood marginal (también llamado evidencia) es la probabilidad de observar los datos $y$ promediado sobre todas los posibles hipótesis (o conjunto de parámetros) $\theta$. Si la oración anterior no es muy clara, no hay problema ya veremos ejemplos que harán más claro este concepto. En general, la evidencia puede ser vista como una simple constante de normalización que en la mayoría de los problemas prácticos puede (y suele) omitirse sin perdida de generalidad. Por lo que el teorema de Bayes suele aparecer escrito como:
\begin{equation} p(\theta \mid y) \propto p(y \mid \theta) p(\theta) \end{equation}El rol de todos estos términos irá quedando más claro a medida que avancemos.
El a posteriori representa todo lo que sabemos de un problema, dado un modelo y un conjunto de datos. Y por lo tanto todas las inferencias estadísticas pueden deducirse a partir de él. Tipicamente esto toma la forma de integrales como la siguiente.
\begin{equation} J = \int \varphi(\theta) \ \ p(\theta \mid y) d\theta \end{equation}Por ejemplo, para calcular la media de $\theta$ deberíamos reemplazar $\varphi(\theta)$, por $\theta$:
\begin{equation} \bar \theta = \int \theta \ \ p(\theta \mid y) d\theta \end{equation}Esto no es más que la definición de un promedio pesado, donde cada valor de $\theta$ es pesado según la probabilidad asignada por el a posteriori.
En la práctica y al resolver problemas computacionales, estas integrales se convierten en simples sumatorias.
El teorema de Bayes es el único estimador usado en estadística Bayesiana. Por lo que conceptualmente la estadística Bayesiana resulta muy simple. Según George Box y Andrew Gelman et al. (2013) la estadística Bayesiana se reduce a tres pasos:
Crear un modelo probabilístico. Los modelos probabilísticos son historias que dan cuenta de como se generan los datos observados (o por observar). Los modelos se expresan usando distribuciones de probabilidad.
Condicionar el modelo a los datos observados a fin de obtener el a posteriori. Usando el teorema de Bayes se actualizan las probabilidades asignadas a priori de acuerdo a los datos observados obteniendose las probabilidades a posteriori.
Criticar el ajuste del modelo generado a los datos y evaluar las consecuencias del modelo. Se puede demostrar que dada la información previa y los datos observados no existe otro mecanismo capaz de generar una mejor inferencia que la estadística Bayesiana. Esto parece maravilloso, pero hay un problema, solo es cierto si se asumen que los datos y el modelo son correctos. En la práctica, los datos pueden contener errores y los modelos son a duras penas aproximaciones de fenómenos reales. Por lo tanto es necesario realizar varias evaluaciones, incluyendo si las predicciones generadas por el modelo se ajustan a los datos observados, si las conclusiones obtenidas tienen sentido dado el marco conceptual en el que uno trabaja, la sensibilidad de los resultados a los detalles del modelo (sobre todo a detalles para los cuales no tenemos demasiada información), etc.
En la práctica la mayoría de los modelos tendrán más de un parámetro, pero al usar software como PyMC3 modelar 1 o 1000 parámetros es más o menos lo mismo. Sin embargo, esos modelos pueden distraernos de los conceptos esenciales, por lo que considero importante comenzar por el caso más sencillo.
Veamos nuestro primer problema:
¿Cuál es la probabilidad de que dicha persona tenga la enfermedad? Antes de seguir leyendo tomate un tiempo y contesta a esta pregunta de forma intuitiva.
A continuación se ve cómo se aplica el teorema de Bayes para dar respuesta a esta pregunta.
$$p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y)}$$$$p(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y \mid \theta)p(\theta) + p(y \mid \theta^c)p(\theta^c)}$$$$p(\theta \mid y) = \frac{0.99 \times 0.001}{0.99 \times 0.001 + 0.05 \times (1−0.001)}$$$$p(\theta \mid y) = \frac{0.00099}{0.00099 + 0.04995}$$$$p(\theta \mid y) \approx 0.02 \approx 2\%$$Si es la primer vez que te encontrás con este problema, las chances son altas de que hayas dado como respuesta un número mucho más alto y cercano a 99%. La mayoría de las personas se dan cuenta que dado que la taza de falsos positivos es del 5% la respuesta tiene que ser menor de 99%, pero fallan en incluir la información a priori, solo 1 de cada mil personas tiene la enfermedad. La importancia del a priori queda claro si uno considera dos casos extremos:
Ésto nos muestra dos cosas, el a priori no puede ser dejado de lado y en general no es buena idea asignar una probabilidad de 1 o 0 a un priori, ya que sin importar lo que nos digan los datos jamás cambiaremos de parecer.
In [2]:
p_apriori_enfermo = 0.001
p_pos_enfermo = 0.99
p_pos_sano = 0.05
(p_pos_enfermo * p_apriori_enfermo) / \
(p_pos_enfermo * p_apriori_enfermo + p_pos_sano * (1-p_apriori_enfermo))
Out[2]:
En el ejemplo anterior el a priori es un número (0.001) en estadística Bayesiana esto equivale a decir que tenemos certeza absoluta sobre el valor del a priori, vale exactamente 0.001 ni un poco más ni un poco menos. Este nivel de precisión es imposible para cualquier problema real, a lo sumo podemos encontrar casos donde el grado de error con el que conocemos una cantidad es muy pequeño en comparación con otros términos de nuestro problema y por lo tanto aproximamos esa cantidad como una constante, aproximación que será válida si no tiene efectos prácticos en nuestras conclusiones.
El ejemplo anterior es un caso sencillo donde la aplicación del teorema de Bayes nos lleva a la respuesta correcta, y donde no hace falta ningún grado de sofisticación matemática para aplicarlo. Pero este ejemplo es apenas un precalentamiento que no hace justicia a la Estadística Bayesiana. Un ejemplo que logra capturar mucho mejor las ideas centrales en estadística Bayesiana es el que veremos a continuación y usaremos el resto de este capítulo y todo el próximo.
A juzgar por la cantidad de ejemplos sobre monedas arrojadas al aires en libros de estadística y probabilidad, pareciera que las monedas son uno de los objetos de estudio centrales de estas disciplinas.
Una de las razones detrás de la ubiquidad de este ejemplo es que las monedas son objetos familiares que facilitan discutir conceptos que de otra forma podrían sonar demasiado abstractos. De todas formas quizá la razón más importante sea que el problema puede ser modelado de forma simple y que muchos problemas reales son conceptualmente similares, de hecho cualquier problema en donde obtengamos resultados binarios (0/1, enfermo/sano, spam/no-spam, etc) puede ser pensado como si estuviéramos hablando de monedas. En definitiva el modelo que veremos a continuación (ejemplificado con monedas) sirve para cualquier situación en la cual los datos observados solo pueden tomar dos valores mutuamente excluyentes. Debido a que estos valores son nominales y son dos, a este modelo se le llama binomial.
En el siguiente ejemplo trataremos de determinar el grado en que una moneda está sesgada. En general cuando se habla de sesgo se hace referencia a la desviación de algún valor (por ejemplo, igual proporción de caras y cecas), pero aquí usaremos el termino sesgo de forma más general. Diremos que el sesgo es un valor en el intervalo [0, 1], siendo 0 para una moneda que siempre cae ceca y 1 para una moneda que siempre cae cara y lo representaremos con la variable $\theta$. A fin de cuantificar $\theta$ arrojaremos una moneda al aire repetidas veces, por practicidad arrojaremos la moneda de forma computacional (¡pero nada nos impide hacerlo manualmente!). Llevaremos registro del resultado en la variable $y$. Siendo $y$ la cantidad de caras obtenidas en un experimento.
Habiendo definido nuestro problema debemos expresarlo en términos del teorema de Bayes,
\begin{equation} p(\theta \mid y) \propto p(y \mid \theta) p(\theta) \end{equation}Donde, como dijimos $\theta = 1$ quiere decir 100% cara y $\theta = 0$ 100% ceca.
Ahora solo restar reemplazar los dos términos a la derecha de la igualdad, el a priori y el likelihood, por distribuciones de probabilidad adecuadas y luego multiplicarlas para obtener el término a la izquierda, el a posteriori. Como es la primera vez que haremos ésto, lo haremos paso a paso y analíticamente. En el próximo capítulo veremos cómo hacerlo computacionalmente.
El a priori lo modelaremos usando una distribución beta, que es una distribución muy usada en estadística Bayesiana. La $pdf$ de esta distribución es:
\begin{equation} p(\theta)= \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, \theta^{\alpha-1}(1-\theta)^{\beta-1} \end{equation}El primer término es una constante de normalización. Por suerte para nuestro problema nos basta con establecer una proporcionalidad, por lo que podemos simplificar esta expresión y escribir la distribución beta de la siguiente forma.
\begin{equation} p(\theta) \propto \theta^{\alpha-1}(1-\theta)^{\beta-1} \end{equation}Hay varias razones para usar una distribución beta para este y otros problemas:
Respecto al último punto, veamos un ejemplo. Supongamos que el experimento de la moneda es realizado por tres personas. Una de ellas dice no saber nada de la moneda por lo tanto a priori todos los valores de $\theta$ son igualmente probables. La segunda persona desconfía de la moneda, ya que sospecha que es una moneda trucada, por lo tanto considera que está sesgada, pero no sabe si hacia cara o hacia ceca. Por último, la tercer persona asegura que lo más probable es que $\theta$ tome un valor alrededor de 0.5 ya que según su experiencia así es como se comportan las monedas. Todas estas situaciones pueden ser modeladas por la distribución beta, como se ve a continuación.
In [3]:
plt.figure(figsize=(10, 3))
x = np.linspace(0, 1, 100)
for ind, (a, b) in enumerate([(1, 1), (0.5, 0.5), (20, 20)]):
y = stats.beta.pdf(x, a, b)
plt.subplot(1, 3, ind+1)
plt.plot(x, y, label='a = %s\nb = %s' % (a, b))
plt.legend(fontsize=12)
plt.tight_layout();
In [4]:
def beta(a, b):
x = np.linspace(0, 1, 130)
beta = stats.beta(a, b)
plt.plot(x, beta.pdf(x))
plt.yticks([])
plt.ylim(0, 6)
interact(beta,
a=ipyw.FloatSlider(min=0.5, max=7, step=0.5, value=2),
b=ipyw.FloatSlider(min=0.5, max=7, step=0.5, value=2));
Habiendo definido el a priori veamos ahora el likelihood. Asumiendo que el resultado obtenido al arrojar una moneda no influye en el resultado de posteriores experimentos (es decir los experimentos son independientes entre sí) es razonable utilizar como likelihood la distribución binomial.
\begin{equation} p(y \mid \theta) = \frac{N!}{y!(N-y)!} \theta^y (1 - \theta)^{N−y} \end{equation}Donde N es la cantidad total de experimentos (monedas arrojadas al aire) e $y$ es la cantidad de caras obtenidas. A los fines prácticos podríamos simplificar la igualdad anterior y convertirla en una proporcionalidad, eliminando el término $\frac{N!}{y!(N-y)!}$ ya que ese término no depende de $\theta$ que es lo que nos interesa averiguar. Por lo que podríamos establecer que:
\begin{equation} p(y \mid \theta) \propto \theta^y (1 - \theta)^{N−y} \end{equation}La elección de esta distribución para modelar nuestro problema es razonable ya que $\theta$ es la chance de obtener una cara al arrojar una moneda y ese hecho ha ocurrido $y$ veces, de la misma forma $1-\theta$ es la chance de obtener ceca lo cual ha sido observado $N-y$ veces.
In [5]:
def binomial(n, θ):
bino = stats.binom(n, θ)
plt.bar(range(n+1), bino.pmf(range(n+1)))
plt.xticks(range(n+1))
plt.ylim(0, 1);
interact(binomial, n=ipyw.IntSlider(min=1, max=10, value=1), θ=ipyw.FloatSlider(min=0, max=1, step=0.05, value=0.5));
Se puede demostrar que siempre que usemos como a priori una función beta y como likelihood una distribución binomial obtendremos como resultado un a posteriori que será nuevamente una distribución beta con los siguientes parámetros:
\begin{equation} p(\theta \mid y) \propto \operatorname{Beta}(\alpha_{a priori} + y, \beta_{a priori} + N - y) \end{equation}Veamos de donde surge este resultado, según el teorema de Bayes el a posteriori es el producto del likelihood y el a priori.
\begin{equation} p(\theta \mid y) \propto p(y \mid \theta) p(\theta) \end{equation}Por lo tanto, en nuestro caso tendremos que:
\begin{equation} p(\theta \mid y) \propto \theta^y (1 - \theta)^{N−y} \theta^{\alpha-1}(1-\theta)^{\beta-1} \end{equation}Reordenando, obtenemos que el a posteriori es:
\begin{equation} p(\theta \mid y) \propto \theta^{\alpha-1+y}(1-\theta)^{\beta-1+N−y} \end{equation}Esto es una distribución Beta (sin considerar la constante de normalización).
Cuando se cumple que para un cierto likelihood la forma funcional del a priori y la del a posteriori coinciden se dice que el a priori es conjugado con el likelihood. Historicamente los problemas en estadística Bayesiana estuvieron restringidos al uso de a prioris conjugados, ya que estos garantizan la tratabilidad matemática del problema, es decir garantizan que es posible obtener una expresión analítica para nuestro problema. En el próximo capítulo veremos técnicas computacionales modernas que permiten obtener a posterioris incluso cuando no se usan a prioris conjugados, lo que ha permitido el resurgimiento de la estadística Bayesiana en las últimas décadas.
Para representar modelos en estadística Bayesiana (y en probabilidad en general) se suele utilizar la siguiente notación
\begin{align} \theta &\sim \operatorname{Beta}(\alpha, \beta) \\ y &\sim \operatorname{Bin}(n=1, p=\theta) \end{align}El símbolo $\sim$ indica que la variable a la izquierda se distribuye según la distribución a la derecha. Entonces podríamos decir que $\mathbf{\theta}$ es una variable aleatoria con distribución beta, y que beta está definida por los parámetros $\alpha$ y $\beta$, este es nuestro a priori. En la siguiente linea tenemos el likelihood el cual está definido por una distribución binomial con parámetros $n=1$ y $p=\theta$. Nótese que tanto $\theta$ como $y$ son vectores, en algunos textos se usa una notación alternativa usando escalares y subíndices como por ejemplo $y_i$.
Gráficamente esto se puede representar usando los diagramas de Kruschke:
En el primer nivel (de arriba hacia abajo) se observa el a priori, luego el likelihood, y por último los datos. Las flechas indican la vinculación entre las partes del modelo y el signo $\sim$ la naturaleza estocástica de las variables.
Bien, ahora que sabemos cómo calcular el a posteriori, lo único que resta es conseguir los datos. En este ejemplo los datos son sintéticos, es decir los obtuve computacionalemnte mediante un generador de números (pseudo)aleatorios, pero bien podrían haber surgido de un experimento con una moneda real.
En el próximo capítulo veremos cómo usar métodos computacionales para computar un a posteriori sin necesidad de derivarlo analíticamente. Esto es lo que haremos para resolver el resto de los problemas del curso. Pero dado que ya nos tomamos el trabajo de derivar analíticamente la expresión para el a posteriori vamos a usar esa expresión. Si miran el código de la siguiente celda verán que la mayoría de las lineas se encargan de dibujar los resultados y no de calcularlos. El cálculo del a posteriori ocurre en las lineas 18, 22 y 26. Cada una de estas lineas computa el a posteriori para cada uno de los a prioris que vimos antes. El cálculo es simple, tan solo se computa el valor del a posteriori (usando la función pdf de la distribución beta provista por SciPy) para 100 puntos igualmente espaciados entre 0 y 1 (linea 7). El loop que empieza en la linea 9 se debe a que exploraremos cómo cambian las distribuciones a posteriori para distinta cantidad de datos (n_intentos). Con una linea negra punteada y vertical se indica el valor real de $\theta$, valor que por supuesto es desconocido en una situación real, pero conocido para mí, ya que los datos son sintéticos.
In [6]:
plt.figure(figsize=(11, 11))
theta_real = 0.35
n_intentos = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
datos = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
beta_params = [(1, 1), (0.5, 0.5), (20, 20)]
dist = stats.beta
x = np.linspace(0, 1, 200)
for idx, N in enumerate(n_intentos):
if idx == 0:
plt.subplot(4, 3, 2)
else:
plt.subplot(4,3, idx+3)
y = datos[idx]
for (a_prior, b_prior), c in zip(beta_params, ('C0', 'C1', 'C2')):
p_theta_dado_y = dist.pdf(x, a_prior + y, b_prior + N - y)
plt.plot(x, p_theta_dado_y, c)
plt.fill_between(x, 0, p_theta_dado_y, color=c, alpha=0.4)
plt.axvline(theta_real, ymax=0.3, color='k')
plt.plot(0, 0, label='%d experimentos\n%d caras' % (N, y), alpha=0)
plt.xlim(0, 1)
plt.ylim(0, 12)
plt.xlabel('θ')
plt.legend()
plt.yticks([])
plt.tight_layout()
La primer figura del panel muestra los a priori, nuestra estimación de $\theta$ dado que no hemos realizado ningún experimento. Las sucesivas nueve figuras muestran las distribuciones a posteriori y se indica la cantidad de experimentos y de caras obtenidas. Además se puede ver una linea negra vertical en 0.35, la cual representa el valor verdadero de $\theta$. Por supuesto que en problemas reales este valor es desconocido.
Este ejemplo es realmente ilustrativo en varios aspectos.
De los ejemplos anteriores debería quedar claro que los a priori influencian los resultados de nuestros cálculos. Esto tiene total sentido si no fuese así no haría falta incluirlos en el análisis y todo sería más simple (aunque nos perderíamos la oportunidad de usar información previa). De los ejemplos anteriores también debería quedar claro que a medida que aumentan los datos (como las tiradas de monedas) los resultados son cada vez menos sensibles al a priori. De hecho, para una cantidad infinita de datos el a priori no tiene ningún efecto. Exactamente cuantos datos son necesarios para que el efecto del a priori sea despreciable varía según el problema y los modelos usados. En el ejemplo de la moneda se puede ver que 50 experimentos bastan para hacer que dos de los resultados sean prácticamente indistinguibles, pero hacen falta más de 150 experimentos para que los 3 resultados se vuelvan practicamente independientes del a priori. Esto es así por que los dos primeros a prioris son relativamente planos, mientras que el tercer a priori concentra casi toda la probabilidad en una región relativamente pequeña. El tercer a priori no solo considera que el valor más probable de $\theta$ es 0.5, si no que considera que la mayoría de los otros valores son muy poco probables. ¿Cómo cambiarían los resultados si hubiéramos usado como a priori $\operatorname{Beta}(\alpha=2, \beta=2)$?
La elección de los a priori puede poner nervioso a quienes se inician en el análisis Bayesiano (o a los detractores de este paradigma). ¡El temor es que los a prioris censuren a los datos y no les permitan hablar por sí mismos! Eso está muy bien, pero el punto es que los datos no saben hablar, con suerte murmuran. Los datos solo tienen sentido a la luz de los modelos (matemáticos y mentales) usados para interpretarlos, y los a prioris son parte de esos modelos.
Hay quienes prefieren usar a priori no-informativos (también conocidos como a priori planos, vagos, o difusos). Estos a priori aportan la menor cantidad posible de información y por lo tanto tienen el menor impacto posible en el análisis. Si bien es posible usarlos, en general hay razones prácticas para no preferirlos. En este curso usaremos a priori ligeramente informativos siguendo las recomendaciones de Gelman, McElreath, Kruschke, y otros. En muchos problemas sabemos al menos algo de los valores posibles que pueden tomar nuestros parámetros, por ejemplo que solo pueden ser positivos, o que están restringidos a sumar 1 o el rango aproximado, etc. En esos casos podemos usar a prioris que introduzcan esta ligera información. En estos casos podemos pensar que la función del a priori es la de mantener las inferencias dentro de límites razonables. Estos a priori se suelen llamar regularizadores.
Por supuesto que también es posible usar a prioris informativos (o fuertes). Hacer esto es razonable solo si contamos con información previa confiable. Esto puede ser ventajoso en casos en que los datos contengan poca información sobre el problema. Si la información no viene por el likelihood (datos), entonces puede venir por el a priori. A modo de ejemplo, en bioinformática estructural es común usar toda la información previa posible (de forma Bayesiana y no-Bayesiana) para resolver problemas. Esto es posible por la existencia de bases de datos que almacenan los resultados de cientos o miles experimentos realizados a lo largo de décadas de esfuerzo (¡No usar esta información sería casi absurdo!). En resumen, si contás con información confiable no hay razón para descartarla, menos si el argumento es algo relacionado con pretender ser objetivo (¡No hay objetividad en negar lo que se sabe!).
Hasta ahora hemos visto que es posible clasificar, aunque sea de forma vaga o aproximada, a los a priori en función de la información que contienen. Pero saber esta clasificación no necesariamente hace las cosas más simples a la hora de elegir un a priori. ¿Acaso no sería mejor eliminar los a prioris de nuestro análisis? Eso haría el asunto mucho mas simple. Bueno, el punto es que desde una perspectiva Bayesiana todos los modelos tienen a prioris, aun cuando no sean explícitos. De hecho muchos resultados de la estadística frecuentista pueden considerarse casos especiales de modelos Bayesianos usando a prioris planos. Volviendo a la figura anterior se puede ver que la moda del a posteriori para la curva azul. Coincide con la estimación (puntual) frecuentista para el valor de $\theta$
\begin{equation} \hat \theta = {{y} \over {N}} \end{equation}Notar que $\hat \theta$ es una estimación puntual (un número) y no una distribución.
Este ejemplo nos muestra que no es posible hacer análisis estadísticos y sacarse los a prioris de encima. Un posible corolario es que es más flexible y transparente especificar los a prioris de forma explícita que esconderlos bajo la cama. Al hacerlo ganamos mayor control sobre nuestro modelo, mayor transparencia y por el mismo precio la estimación de la incertidumbre con la que se estima cada parámetro.
Por último, hay que recordar que el modelado estadístico (como otras formas de modelado) es un proceso iterativo e interactivo. Nada nos impide usar más de un a priori (o un likelihood) si así lo quisiéramos. Una parte importante del modelado es la de cuestionar los supuestos y los a prioris son simplemente un tipo de supuestos (como lo son los likelihoods). Si tuvieramos más de un a priori razonable podríamos realizar un análisis de sensibilidad, es decir evaluar como cambian los resultados con los a prioris, podríamos llegar a la conclusión que para un rango amplio de a prioris ¡los resultados no varían! Más adelante veremos varias herramientas para comparar distintos modelos.
Dado que los a prioris tienen un papel central en la estadística Bayesiana, seguiremos discutiéndolos a medida que vayamos viendo problemas concretos. Por lo que si esta discusión no ha aclarado todas tus dudas y seguís algo confundido, mejor mantener la calma y no preocuparse demasiado, este tema ha sido motivo de discusión y confusión durante décadas ¡y la discusión todavía continua!
En general la distribución más familiar para la mayoría de las personas es la distribución Gaussiana, como esta distribución está definida por dos parámetros, la media y la dispersión de ese valor medio, suele resultarnos natural pensar las distribuciones en esos términos. Si queremos expresar la distribución beta en función de la media y la dispersión podemos hacerlo de la siguiente forma:
\begin{align} \alpha &= \mu \kappa \\ \beta &= (1 − \mu) \kappa \end{align}donde $\mu$ es la media y $\kappa$ es un parámetro llamado concentración. Por ejemplo si $\mu=0.5$ y $\kappa=40$, tenemos que:
\begin{align} \alpha &= 0.5 \times 40 &= 20 \\ \beta &= (1-0.5) \times 40 &= 20 \end{align}$\kappa$ se puede interpretar como la cantidad de experimentos si/no que realizamos dándonos como resultado la media $\mu$. Es decir el a priori no sesgado (verde) equivale a haber arrojado una moneda 40 veces y haber obtenido como media 0.5. Es decir que si usamos ese a priori recién al observar 40 experimentos si/no, los datos tendrán el mismo peso relativo que el a priori, por debajo de este número el a priori contribuye más que los datos al resultado final y por encima menos. El a priori azul (uniforme) equivale a haber observado a la moneda caer una vez cara y otra vez ceca ($\kappa = 2$). Cuando $\kappa < 2$, la cosa se pone un poco extraña, por ejemplo el a priori sesgado (anaranjado) equivale a haber observado una sola moneda ($\kappa = 1$) pero en una especie de (a falta de mejor analogía) ¡superposición cuántica de estados!
El resultado de un análisis Bayesiano es siempre una distribución de probabilidad. En el caso de la moneda esto es evidente, y en el caso del diagnostico es menos claro ya que la distribución es discreta y solo puede tomar dos valores.
A la hora de comunicar los resultados de un análisis Bayesiano, lo más informativo es reportar la distribución completa, aunque esto no siempre es posible o deseable, por ejemplo el a posteriori de una distribución multidimensional es imposible de dibujar en papel. En general, se suele recurrir a distintas medidas que resumen el a priori, por ejemplo reportando la media de la distribución a posteriori. Algo un poco más informativo es reportar además un intervalo de credibilidad. Existen varios criterios para definir intervalos de credibilidad, el que usaremos en este curso (y que también es ampliamente usado en la literatura) es lo que se conoce como intervalo de más alta densidad y nos referiremos a él por su sigla en ingles, HPD (Highest Posterior Density interval). Un HPD es el intervalo, más corto, que contiene una porción fija de la densidad de probabilidad, generalmente el 95% (aunque otros valores como 90% o 50% son comunes). Cualquier punto dentro de este intervalo tiene mayor densidad que cualquier punto fuera del intervalo. Para una distribución unimodal, el HPD 95 es simplemente el intervalo entre los percentiles 2,5 y 97,5.
ArviZ es un paquete de Python para análisis exploratorio de modelos Bayesianos. ArviZ provee de funciones que facilitan el resumir el a posteriori. Por ejemplo `plot
El resultado de un análisis Bayesiano es siempre una distribución de probabilidad. En el caso de la moneda esto es evidente, y en el caso del diagnostico es menos claro ya que la distribución es discreta y solo puede tomar dos valores.
A la hora de comunicar los resultados de un análisis Bayesiano, lo más informativo es reportar la distribución completa, aunque esto no siempre es posible o deseable, por ejemplo el a posteriori de una distribución multidimensional es imposible de dibujar en papel. En general, se suele recurrir a distintas medidas que resumen el a priori, por ejemplo reportando la media de la distribución a posteriori. Algo un poco más informativo es reportar además un intervalo de credibilidad. Existen varios criterios para definir intervalos de credibilidad, el que usaremos en este curso (y que también es ampliamente usado en la literatura) es lo que se conoce como intervalo de más alta densidad y nos referiremos a él por su sigla en ingles, HPD (Highest Posterior Density interval). Un HPD es el intervalo, más corto, que contiene una porción fija de la densidad de probabilidad, generalmente el 95% (aunque otros valores como 90% o 50% son comunes). Cualquier punto dentro de este intervalo tiene mayor densidad que cualquier punto fuera del intervalo. Para una distribución unimodal, el HPD 95 es simplemente el intervalo entre los percentiles 2,5 y 97,5.
ArviZ es un paquete de Python para análisis exploratorio de modelos Bayesianos. ArviZ provee de funciones que facilitan el resumir el a posteriori. Por ejemplo plot_posterior puede ser usado para generar un gráfico con la media y HPD. En el siguiente ejemplo en vez de un a posteriori de un ejemplo real estamos usando datos generados al azar según una distribución beta.
In [7]:
np.random.seed(1)
mock_posterior = stats.beta.rvs(5, 11, size=1000) # az.convert_to_inference_data({'θ':stats.beta.rvs(5, 11, size=1000)}, chains=1)
az.plot_posterior(mock_posterior, figsize=(8, 4), textsize=15);
Ahora que estamos aprendiendo que es un HPD por primera vez y antes de que automaticemos el concepto conviene aclarar un par de puntos.
La elección automática de 95% (o cualquier otro valor) es totalmente arbitraria. En principio no hay ninguna razón para pensar que describir el a posteriori con un HPD 95 sea mejor que describirlo con un HPD 98 o que no podamos usar valores como 87% o 66%. El valor de 95% es tan solo un accidente histórico. Como un sutil recordatorio de esto ArviZ usa por defecto el intervalo de 94%.
Un intervalo de credibilidad (que es Bayesiano) no es lo mismo que un intervalo de confianza (que es frecuentista). Un intervalo de confianza es un intervalo que se define según un nivel de confianza, en general del 95%. Un intervalo de confianza se construye de tal forma que si repitiéramos infinitas veces un experimento obtendríamos que la proporción de intervalos que contienen el valor verdadero del parámetro que nos interesa coincide con el nivel de confianza estipulado. Contra-intuitivamente esto no es lo mismo que decir que un intervalo en particular tiene una probabilidad $x$ de contener el parámetro (esto sería la definición de un intervalo de credibilidad, que es Bayesiano). De hecho, un intervalo de confianza en particular contiene o no contiene al valor, la teoría frecuentista no nos deja hablar de probabilidades de los parámetros, ya que estos tienen valores fijos. Si no queda clara la diferencia no te hagas problema, la diferencia entre estos dos conceptos suele ser tan difícil de entender que en la práctica estudiantes y científicos por igual interpretan los intervalos de confianza (frecuentistas) como intervalos de credibilidad (Bayesianos).
Si bien desde la perspectiva Bayesiana podemos afirmar que un intervalo de credibilidad nos permite asegurar que la probabilidad de un parámetro está acotado en cierto rango. Siempre hay que tener presente que dicha afirmación es correcta SOLO en sentido teórico. Es decir, solo si todos los supuestos contenidos en el modelo son ciertos. Una inferencia es siempre dependiente de los datos y modelos usados.
Las pruebas predictivas a posteriori son una forma de evaluar el modelo. Una vez obtenido el modelo Bayesiano se usa el a posteriori para generar datos $\tilde{y}$, es decir datos predichos condicionados por los valores estimados de $\theta$ y por los datos ya observados $y$.
\begin{equation} p(\tilde{y} \mid y) = \int p(\tilde{y} \mid \theta) p(\theta \mid y) d\theta \end{equation}Los datos generados son predictivos ya que son los datos que se esperaría ver por ejemplo en un futuro experimento, es decir son variables no observadas pero potencialmente observables. La prueba consiste en comparar los datos observados con los datos predichos a partir del a posteriori.
Las pruebas predictivas a posteriori son pruebas de auto-consistencia. Este ejercicio nos permite evaluar si el modelo es razonable, la idea general no es determinar si un modelo es correcto o no ya que como dijo George Box "todos los modelos están equivocados, pero algunos son útiles". El grado de confianza en la verosimilitud de los modelos ciertamente es distinta entre practicantes de distintas disciplinas científicas, en disciplinas como física cuando se estudian sistemas relativamente simples bajo condiciones experimentales extremadamente controladas y haciendo uso de teorías fuertes, es probable que se le asigne un alto grado de confianza a ciertos modelos. Pero esto no suele ser cierto en disciplinas como ciencias sociales o biología (aunque sospecho que la variabilidad encontrada en biología ¡es muy alta!). En el caso de contar a prioris muy informativos la evaluación de un modelo también puede ser usado para evaluar si los propios datos son razonables, indicando que tal vez sea necesario conseguir nuevos datos o revisar como se obtuvieron los datos o como se procesaron.
En definitiva la principal utilidad de las pruebas predictivas a posteriori debería ser el permitirnos dar una segunda mirada, crítica, al modelo y tratar de entender la razón de discrepancias sistemáticas (si las hubiera), estas discrepancias nos pueden llevar a entender mejor los límites del modelo, abandonar el modelo por completo o tal vez mejorarlo.
Si bien se han desarrollado métodos formales o cuantitativos para realizar pruebas predictivas a posteriori, una aproximación que suele ser más informativa y simple de interpretar es realizar gráficas, que es lo que iremos viendo en los próximos capítulos.
Empezamos este capítulo con una breve discusión sobre el modelado estadístico y la teoría de la probabilidad y teorema de Bayes que se deriva de ella. Luego utilizamos el problema de la moneda como una excusa para introducir aspectos básicos del modelado Bayesiano y el análisis de datos. Utilizamos este ejemplo clásico para transmitir algunas de las ideas más importantes de las estadística bayesiana, fundamentalmente el uso de distribuciones de probabilidad para construir modelos y representar la incertidumbre. Tratamos de desmitificar el uso de los a prioris dandoles el mismo estatus epistemológico-metodológico que otros elementos que forman parte del proceso de modelado e inferencia, como el likelihood o incluso meta-preguntas, como por qué estamos tratando de resolver un problema en particular. Concluimos el capítulo con una breve descripción de cómo interpretar y comunicar los resultados de un análisis bayesiano.
La siguiente figura, inspirada en una figura de Sumio Watanabe resume el flujo de trabajo Bayesiano tal cual se describió en este capítulo.
Es importance notar que mientras la distribución a posteriori es una distribution sobre los parámetros en un modelo, la distribución predictiva a posteriori es una distribución sobre los datos (predichos).
De las siguientes expresiones cual(es) se corresponde(n) con el enunciado "la probabilidad de lluvia dado que es 25 de Mayo de 1810"?
Enuncie con palabras cada una de las expresiones del punto anterior.
In [8]:
def posterior_grid(grilla=10, a=1, b=1, caras=6, tiradas=9):
grid = np.linspace(0, 1, grilla)
prior = stats.beta(a, b).pdf(grid)
likelihood = stats.binom.pmf(caras, tiradas, grid)
posterior = likelihood * prior
posterior /= posterior.sum()
_, ax = plt.subplots(1, 3, sharex=True, figsize=(16, 4))
ax[0].set_title('caras = {}\ntiradas = {}'.format(caras, tiradas))
for i, (e, e_n) in enumerate(zip([prior, likelihood, posterior], ['prior', 'likelihood', 'posterior'])):
ax[i].set_yticks([])
ax[i].plot(grid, e, 'o-', label=e_n)
ax[i].legend(fontsize=14)
interact(posterior_grid, grilla=ipyw.IntSlider(min=2, max=100, step=1, value=15), a=ipyw.FloatSlider(min=1, max=7, step=1, value=1), b=ipyw.FloatSlider(
min=1, max=7, step=1, value=1), caras=ipyw.IntSlider(min=0, max=20, step=1, value=6), tiradas=ipyw.IntSlider(min=0, max=20, step=1, value=9));