Esta notebook fue creada originalmente como un blog post por Raúl E. López Briega en Matemáticas, análisis de datos y python. El contenido esta bajo la licencia BSD.
La estadística es la ciencia que relaciona los datos con las preguntas que nos interesa responder. Esto incluye elaborar métodos para recopilar los datos pertinentes a cada pregunta, métodos para resumir y mostrar los datos de forma tal que arrojen luz sobre las preguntas a responder; y los métodos que nos permitan dibujar las respuestas a cada pregunta y que sean apoyadas por los datos. Cuando trabajamos con datos, debemos tener en cuenta que casi siempre contienen incertidumbre. Esta incertidumbre puede surgir de la selección de los elementos que tomamos para nuestras mediciones, o puede surgir de la variabilidad misma del proceso de medición. En cierto sentido, podemos decir que obtener conclusiones generales de los datos es la base para aumentar el conocimiento que tenemos sobre el mundo,y es la base de toda investigación científica.
El método científico busca relaciones de causa y efecto entre una variable experimental y una variable de resultado. Es decir, cómo cambiando el valor de la variable experimental resulta en un cambio en la variable de resultado. En muchas oportunidades este proceso puede ser tan simple o complejo como desarrollar modelos matemáticos para estas relaciones; pero siempre es necesario aislar al experimento de factores externos que podrían afectar los resultados experimentales. Estos factores externos deben ser controlados de forma tal que el verdadero efecto de la variable experimental en la variable de resultado pueda ser determinado con cierto grado de seguridad. En ciencias como biología, medicina, ingeniería, tecnología y ciencias sociales no es tan fácil identificar estos factores relevantes que deben ser controlados. En esos campos se necesita una forma diferente de controlar a los factores externos; es aquí dónde la inferencia estadística se vuelve una herramienta fundamental para apoyar a la investigación científica.
Los métodos estadísticos de inferencia pueden utilizarse cuando hay una variabilidad aleatoria en los datos. Esto puede extender el método científico a situaciones en las que los factores externos relevantes no pueden ser identificados. Como no podemos identificar estos factores externos, no podemos controlarlos directamente. La falta de control directo significa que los factores externos afectarán los datos, por lo que existirá el peligro de que podamos extraer conclusiones erróneas del experimento debido a estos factores externos que no se pueden controlar. La estadística nos ofrece la importante idea de la asignación aleatoria o al azar para hacer frente a esta posibilidad. Los factores externos no identificados pueden ser promediados asignando aleatoriamente cada unidad a un grupo de tratamiento, o a un grupo de control. La asignación aleatoria nos permite controlar estadísticamente los factores externos por medio de promediar sus efectos; aunque esto no evitará que las conclusiones estadísticas tengan un cierto grado de incertidumbre o error debido a esta variabilidad. La asignación aleatoria no sólo nos ayuda a reducir la incertidumbre debido a los factores externos, sino que también nos permite medir la cantidad de incertidumbre que permanece en las conclusiones del experimento apoyados por un modelo de probabilidad.
Dentro de este contexto, podemos encontrar a la inferencia bayesiana, la cual se caracteriza justamente por tratar a los elementos desconocidos en forma aleatoria; y nos permite ir adaptando el modelo a medida que nos vamos encontrando con nueva evidencia. Es decir, que vamos a comenzar con una distribución previa de los datos, la cual puede ser elegida al azar o en base a nuestra experiencia; y luego vamos a ir adaptando el modelo a medida que vamos agregando nueva evidencia hasta alcanzar una distribución posterior. Esta distribución posterior será la que mejor se adapte a la información con que contábamos al principio, mas la nueva evidencia que hemos adquirido y va a representar el grado de credibilidad que le asignamos a nuestra hipótesis. La importancia de la inferencia bayesiana la vamos a encontrar sobre todo cuando no contamos con una gran cantidad de información para comprobar nuestro modelo; nos va a permitir expresar la incertidumbre sobre el resultado de lo que estamos modelando. Pero para entender en mayor profundidad a la inferencia bayesiana comencemos por el principio.
Thomas Bayes fue un ministro presbiteriano y matemático inglés que estudió la relación íntima que existe entre la probabilidad, la predicción y el progreso científico. Su trabajo se centró principalmente en cómo formulamos nuestras creencias probabilísticas sobre el mundo que nos rodea cuando nos encontramos con nuevos datos o evidencias. El argumento de Bayes no es que el mundo es intrínsecamente probabilístico o incierto, ya que él era un creyente en la divina perfección; sino que aprendemos sobre el mundo a través de la aproximación, acercándonos cada vez más a la verdad a medida que recogemos más evidencias. Este argumento lo expresó matemáticamente a través de su famoso teorema:
$$P(H|D) = \frac{P(D|H)P(H)}{P(D)} $$En donde:
Si los fundamentos filosóficos del teorema de Bayes son sorprendentemente ricos, sus matemáticas son increíblemente simples. En su forma más básica, no es más que una expresión algebraica con tres variables conocidas y una incógnita; y que trabaja con probabilidades condicionales, nos dice la probabilidad de que una hipótesis $H$ sea verdadera si algún evento $D$ ha sucedido. El teorema de Bayes es útil porque lo que normalmente sabemos es la probabilidad de los efectos dados las causas, pero lo que queremos saber es la probabilidad de las causas dadas los efectos. Por ejemplo, sabemos qué porcentaje de pacientes con gripe tiene fiebre, pero lo que realmente queremos saber es la probabilidad de que un paciente con fiebre tenga gripe. El teorema de Bayes nos permite ir de uno a otro con suma facilidad.
Para que quede más claro, ilustremos el teorema con un simple ejemplo del diagnostico médico, una de las aplicaciones dónde más éxito ha tenido la inferencia bayesiana. Supongamos que nos hicimos un estudio y nos ha dado positivo para una rara enfermedad que solo el 0.3 % de la población tiene; la tasa de efectividad de este estudio es del 99 %, es decir, que solo da falsos positivos en el 1 % de los casos. ¿Cuán probable es que realmente tengamos la enfermedad?. En un principio, nos veríamos tentados a responder que hay un 99 % de probabilidad de que tengamos la enfermedad; pero en este caso nos estaríamos olvidando del concepto importante del a priori, sabemos con anterioridad que la enfermedad es extremadamente rara (solo el 0.3 % la tiene), si incluimos esta información previa a nuestro cálculo de probabilidad y aplicamos el teorema de Bayes podemos llegar a una conclusión totalmente distinta.
$$ P(\text{ rara enfermedad | positivo}) = \frac{P(\text{ positivo | rara enfermedad})P( \text{rara enfermedad})}{P(\text{positivo})}$$
In [1]:
# Ejemplo simple teorema de Bayes aplicado a estimación de un sólo parámetro.
a_priori = 0.003
likelihood = 0.99
evidencia = 0.01
a_posteriori = likelihood * a_priori / evidencia
a_posteriori
Out[1]:
Como vemos, luego de aplicar el teorema de Bayes llegamos a la conclusión de que en realidad nuestra probabilidad de estar realmente enfermo es de sólo 30 % y no de 99 %, ya que podemos ser uno de los falsos positivos del estudio y la enfermedad es realmente muy rara. Como este ejemplo demuestra, la inclusión del a priori es sumamente importante para la inferencia bayesiana, por lo cual también debemos ser sumamente cuidadosos a la hora de elegirla. Cuando nuestra a priori es fuerte, puede ser sorprendentemente resistente frente a nuevas evidencias.
Los problemas de monedas son clásicos cuando hablamos de probabilidad y estadística, nos permiten ejemplificar conceptos abstractos de forma simple. Asimismo, pueden ser muchas veces conceptualmente similares a situaciones reales, de hecho cualquier problema en donde obtengamos resultados binarios, 0/1, enfermo/sano, spam/no-spam, puede ser pensado como si estuviéramos hablando de monedas.
In [2]:
# <!-- collapse=True -->
# importando modulos necesarios
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import seaborn as sns
import pymc3 as pm
import theano.tensor as tt
np.random.seed(1984) #replicar random
%matplotlib inline
In [3]:
n_experimentos = 4
theta_real = .35 # en una situación real este valor es desconocido
datos = stats.bernoulli.rvs(theta_real, size=n_experimentos)
datos
Out[3]:
In [4]:
with pm.Model() as nuestro_primer_modelo:
# a priori
theta = pm.Beta('theta', alpha=1, beta=1)
# likelihood
y = pm.Bernoulli('y', p=theta, observed=datos)
#y = pm.Binomial('theta',n=n_experimentos, p=theta, observed=sum(datos))
In [5]:
with nuestro_primer_modelo:
trace = pm.sample(1000)
In [6]:
burnin = 0 # por ahora lo fijamos en 0, es decir nada de burnin
pm.traceplot(trace[burnin:], lines={'theta':theta_real});
# lines={'theta':theta_real} permite mostrar el valor de theta_real (el valor "correcto"),
# como una linea roja
In [7]:
pm.df_summary(trace)
Out[7]:
In [8]:
import numpy as np
import matplotlib.pyplot as plt
# Initialize random number generator
np.random.seed(123)
# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]
# Size of dataset
size = 100
# Predictor variable
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2
# Simulate outcome variable
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma
In [9]:
%matplotlib inline
fig, axes = plt.subplots(1, 2, sharex=True, figsize=(10,4))
axes[0].scatter(X1, Y)
axes[1].scatter(X2, Y)
axes[0].set_ylabel('Y'); axes[0].set_xlabel('X1'); axes[1].set_xlabel('X2');
In [10]:
from pymc3 import Model, Normal, HalfNormal
basic_model = Model()
with basic_model:
# Priors for unknown model parameters
alpha = Normal('alpha', mu=0, sd=10)
beta = Normal('beta', mu=0, sd=10, shape=2)
sigma = HalfNormal('sigma', sd=1)
# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2
# Likelihood (sampling distribution) of observations
Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
In [4]:
from pymc3 import find_MAP
map_estimate = find_MAP(model=basic_model)
print(map_estimate)
In [5]:
from pymc3 import NUTS, sample
from scipy import optimize
with basic_model:
# draw 500 posterior samples
trace = sample(500)
In [6]:
trace['alpha'][-5:]
Out[6]:
In [7]:
from pymc3 import Slice
with basic_model:
# obtain starting values via MAP
start = find_MAP(fmin=optimize.fmin_powell)
# instantiate sampler
step = Slice(vars=[sigma])
# draw 5000 posterior samples
trace = sample(5000, step=step, start=start)
In [8]:
from pymc3 import traceplot
traceplot(trace);
In [9]:
from pymc3 import summary
summary(trace)
In [ ]:
In [ ]: