In [1]:
%matplotlib inline
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from ipywidgets import interact
plt.style.use('seaborn-darkgrid')
En muchos casos describir los datos y graficarlos puede ser suficiente para nuestros propósitos. Pero en otros casos esto es solo el primer paso de un análisis que intenta no solo describir los datos si no comprender los procesos o mecanismos que dieron origen a las observaciones. Esto se conoce como estadística inferencial y es el tema central de este capítulo.
La estadística es una forma de modelado. Los modelos son una descripción, simplificada, de un sistema bajo estudio y no solo una descripción de los datos. Los modelos se construyen intentando capturar de la mejor forma posible la información relevante, por lo que un modelo más complejo no necesariamente es mejor que uno simple. Decenas, cientos o miles de datos pueden pasar a ser descriptos por unos pocos parámetros que definen al modelo, por lo que un modelo es además una forma de compresión de información.
Dentro de la estadística Inferencial se suelen distinguir dos grande paradigmas:
En los últimos años se ha desarrollado, además, un conjunto de técnicas agrupadas bajo el rótulo de machine learning. El machine learning es algo así como una reinvención/re-descubrimiento de la estadística por parte de informáticos más preocupados en resolver problemas que en demostrar teoremas.
Una de las diferencias, algo caricaturizada (y exagerada), entre la estadística inferencial y el machine learning es el énfasis en el modelo. En machine learning se usan modelos que en general son de caja negra. Es decir, los parámetros del modelo no son de interés principal, lo importante es predecir. En cambio la estadística inferencial pone énfasis en entender los parámetros del modelo como medio para entender el problema y eventualmente (aunque no siempre) predecir el comportamiento del sistema. Esta linea que separa ambas aproximaciones se está tornando cada vez más difusa, el flujo de ideas entre estas dos disciplinas (que por mucho tiempo se desarrollaron en paralelo), es cada vez mayor. Los métodos de machine learning usan muchas ideas de la estadística. Además, es cada vez más grande el interés por parte de la comunidad de machine learning en desarrollar métodos más fáciles de interpretar y también de incorporar conceptos de la estadística Bayesiana. En el fondo estas tres grandes ramas de la inferencia tienen un mismo sustento teórico que es la teoría de probabilidad.
Un modelo mental muy útil en estadística consiste en considerar que existe una población que es es el conjunto de objetos/variables que me interesan. En general la población no es algo a lo que tengamos acceso directo, con suerte solo tenemos acceso a un subconjunto de esa población al que llamamos muestra (otros nombres posibles son datos, dataset, observaciones, "data").
El proceso de inferencia consiste en poder estimar propiedades de la población a partir de la muestra. Es común hablar de dos grandes métodos de inferencia:
Al analizar datos, se puede pensar en dos grandes alternativas:
Ambas aproximaciones pueden encararse desde una posición frecuentista o Bayesiana. En este capítulo veremos la forma frecuentista y dejaremos la aproximación Bayesiana para el tercer curso de esta especialidad.
Dentro del paradigma frecuentista la primera de estas estrategias es la que tradicionalmente se ha enseñado y divulgado con más énfasis. Mientras que la segunda viene ganando adeptos y se la ha llamado "nueva estadística" (aunque las ideas distan mucho de ser realmente novedosas).
Las pruebas de hipótesis no suelen ser demasiado comunes en Ciencia de Datos, pero de todas formas merece la pena visitarlas, aunque sea brevemente, ya que suelen ser objeto de muchas confusiones y malentendidos.
La prueba de hipótesis,también llamada test de hipótesis o prueba de significación, es un procedimiento para establecer si una propiedad dada de una población es compatible con lo observado en una muestra de dicha población.
En la práctica, como la población es desconocida, lo que se hace es pensar una hipótesis y partir de la misma generar una distribución de muestreo. Sobre esta distribución se calculan probabilidades. Luego se compara lo observado con lo calculado y esa diferencia se expresa en términos de la significancia estadística.
Este proceso puede hacerse tanto analíticamente como de forma númerica/computacional. Si alguna vez calcularon algo en estadística que luego tenían que buscar en una tabla llena de números, entonces ustedes seguro que hicieron la versión analítica. Nosotros veremos ambas formas.
De forma general y resumida, tenemos:
Definimos un estadístico $T_{\mathcal{D}}$ es decir una medida cuantitativa que describe nuestros datos $\mathcal{D}$
Imaginamos una hipótesis nula $H_{0}$, es decir la hipótesis del no-efecto (la no-diferencia) y asumiendo $H_{0}$ calculamos $T_{\mathcal{\tilde D}}$
Calculamos $\text{valor p} \triangleq p(T_{\mathcal{\tilde D}} \ge T_{\mathcal{D}} ) \mid H_{0})$. Es decir la probabilidad de obtener un estadístico al menos tan extremo como el observado, asumiendo $H_{0}$ como cierto.
Evaluamos si el $\text{valor p}$ es pequeño o grande. Para esto usamos un valor predefinido, el cual suele ser $\alpha = 0.05$. Si $p < \alpha$ entonces se dice que "rechazamos $H_0$" si en cambio $p \ge \alpha$ se dice que "fallamos en rechazar $H_0$".
La lógica de este procedimiento es similar a la prueba por contradicción usada en lógica y matemática. La cual esencialmente sigue el siguiente razonamiento: para probar X, asumimos por un momento que X es falso, si ese supuesto nos guía a un resultado contradictorio entonces concluimos que X debe ser cierto.
De forma similar (aunque no idéntica) en una prueba de hipótesis, para evaluar si el efecto que observamos/medimos existe, asumimos una hipótesis nula que dice que tal efecto no existe, y basados en ese supuesto calculamos la probabilidad de observar nuestros datos. Ese es el $\text{valor p}$. Si el $\text{valor p}$ es bajo (según un criterio predefinido) concluimos que los datos no parecen estar de acuerdo con la hipótesis nula. Es importante destacar que al usar este procedimiento NO podemos establecer que la hipótesis nula sea cierta o falsa, lo que estamos calculando es si los datos parecen tener sentido cuando asumimos la hipótesis nula como cierta.
Para poner estas ideas en un plano más concreto veamos un ejemplo tan trillado como útil. Supongamos que arrojamos al aire 200 veces una misma moneda y observamos 90 caras y 110 cecas, y entonces nos preguntamos ¿Es compatible una moneda que cae la mitad de las veces cara y la mitad ceca con estos datos?
El primer paso es calcular un estadístico, $T_{\mathcal{D}}$, que para este ejemplo bien podría ser la cantidad de caras $z$ dividido la cantidad de tiradas $n$. Por lo tanto nuestro estadístico está restringido a valores entre 0 y 1.
El segundo paso es imaginar una hipótesis nula. Dada la pregunta suena razonable usr como hipótesis nula tenemos el mismo número de caras que de cecas, es decir $H_{0}$ implica que $z$ es la mitad de $n$ y por lo tanto $T_{\mathcal{\tilde D}} = 0.5$.
El tercer paso es calcular el $\text{valor p}$. Para ellos vamos a usar Python, vamos a crear una lista monedas_obs que contiene 200 elementos, 90 1 (representando las caras) y 110 0 (representando las cecas).
In [2]:
monedas_obs = [1] * 90 + [0] * 110
#monedas_obs = [1] * 32 + [0] * 44
# una forma alternativa a la linea de arriba, pero usando NumPy
# monedas_obs = np.repeat((1, 0), (90, 110))
n = len(monedas_obs)
z = sum(monedas_obs)
Ahora que tenemos los datos vamos a calcular el $\text{valor p}$ mediante una simulación. Para ello
In [3]:
dist_muestreo = []
for _ in range(10000):
monedas_sim = np.random.choice([0, 1], size=n)
z_sim = np.sum(monedas_sim)
T = z_sim / n
dist_muestreo.append(T)
y = np.array(dist_muestreo)
T_obs = z / n
valor_p = (np.sum([y <= T_obs]) + np.sum([y >= 1-T_obs])) / len(y)
kde = stats.gaussian_kde(y)
x = np.linspace(y.min(), y.max(), 150)
density = kde.pdf(x)
plt.plot(x, density)
plt.fill_between(x[x <= T_obs], density[x <= T_obs], color='C1',
label='p = {:.3f}'.format(valor_p))
plt.fill_between(x[x >= 1-T_obs], density[x >= 1-T_obs], color='C1')
plt.yticks([])
plt.legend()
plt.xlabel(r'$T_{\mathcal{\tilde D}}$', fontsize=14);
Recapitulemos; lo que acabamos de hacer es calcular la probabilidad de obtener fracciones de caras respecto del total de las tiradas (curva azul), este número va entre 0 (ninguna cara observada en las 200 tiradas) y 1 (100% de caras en las 200 tiradas). La curva azul se llama distribución de muestreo (sampling distribution). Y es la distribución de nuestro estadístico $T_{\mathcal{\tilde D}}$ para una muestra finita de $n=200$ asumiendo $H_{0}$.
Esta distribución es central en estadística frecuentista. De hecho a partir de ella hemos calculado el $\text{valor p}$, que se corresponde con la suma de las dos áreas naranja en las colas de la curva azul.
¿Por qué tenemos dos colas?
Esto deriva de la pregunta que estamos intentando contestar: "Son nuestros datos compatibles con una moneda no sesgada, es decir con una moneda que cae la mitad de las veces cara y la mitad ceca". Dado que es posible obtener una moneda sesgada "para un lado (las caras) o para el otro (las cecas)" el $\text{valor p}$ tiene en cuenta ambos casos. Esto se llama prueba de dos colas, si hubiéramos argumentado que solo nos interesa un caso hubieramos calculado uno solo de esos casos y tendríamos una prueba de una cola. Ambos análisis son totalmente válidos, cual usar depende de la pregunta a contestar.
¿Cómo se interpreta el $\text{valor p}$?
Bueno, desde el punto de vista frecuentista las probabilidades viven en asintopia, un mundo donde las probabilidades son frecuencias obtenidas en el límite de infinitas repeticiones de un mismo experimento. Entonces la interpretación sería:
Si arrojáramos 200 veces al aire una moneda NO SESGADA y esto lo repetimos infinitas veces entonces el ~9% de las veces obtendremos un sesgo al menos tan grande como el observado para caras. Y el ~18% obtendremos un resultado al menos tan sesgado como el observado tanto para caras como para cecas.
Es decir hemos encontrado la respuesta al pregunta ¿Cuán esperable es este resultado para una moneda que cae la mitad de las veces cara?
Otra pregunta natural sería ¿Es realmente diferente este resultado de $H_{0}$? Para dar respuesta a esta pregunta se suele usar el concepto (o "confupto", es decir concepto muy confuso) de significancia estadística, que consiste, como ya adelantamos, en comparar el $\text{valor p}$ con un valor predeterminado, denominado usualmente $\alpha$. Se ha vuelto común utilizar ciertos valores de $\alpha$ como referencia cuasi-universales y tablas de correspondencia ente $\text{valor p}$ y niveles de significancia, por ejemplo:
En este caso como p > 0.05 diríamos: "Se falla en rechazar la hipótesis nula con un $p \simeq 0.18$ y un nivel de significancia de 0.05".
¿De donde salen estos valores predeterminados?
La verdad, de una (con)fusión de dos métodos "rivales" de inferencia, junto con un accidente histórico que dejo cuasi-congelado el valor de $\alpha$ en 0.05.
Fuera de este detalle histórico, la intención de fijar un valor de $\alpha = 0.05$ (o cualquier otro) es la de controlar el nivel de errores de tipo I. Es decir las falsas alarmas o más formalmente controlar la cantidad de veces que estamos dispuestos a rechazar una hipótesis nula verdadera (es decir una que NO deberíamos haber rechazado). Lamentablemente se ha vuelto muy común que estos límites se utilicen para tomar decisiones automáticas por ej si un resultado deber ser publicado. El límite entonces lo fija un criterio editorial (aunque auto-impuesto por los propios científicos) que dependerá de la revista/disciplina, siendo 0.05 común por ej en psicología y en varias ramas de la biología. Esto es claramente absurdo ya que ese límite debería depender del problema que se está estudiando y de detalles que los especialistas en el tema entienden y que en caso de ser controvertidos deberían quedar explícitos para su discusión.
Un par de puntos a tener en cuenta:
Bajo el paradigma frecuentista no es posible aceptar una hipótesis nula. Por eso los frecuentistas dicen "fallamos en rechazar $H_{0}$".
El computo del $\text{valor p}$ asume que $H_{0}$ es cierta. Por lo tanto el $\text{valor p}$ NO es la probabilidad que $H_{0}$ sea cierta, tampoco es la probabilidad de $no H_{0}$ (la alternativa) sea cierta, ni siquiera es $p(\mathcal{D} \mid H_{0})$.
Bajo el paradigma frecuentista no es posible preguntarse sobre la probabilidad que un parámetro tome tal o cual valor. Los parámetros tienen valores fijos (aunque desconocidos). Por lo tanto NO podemos averiguar $p(z = 0.5 \mid \mathcal{D})$, es decir ¿Cual es la probabilidad, dado los datos, que nuestra moneda NO esté sesgada?
Los $\text{valores p}$, para este problema, pueden ser calculados de forma analítica, basta aplicar la distribución binomial.
¿Por qué esta distribución?
Precisamente por que modela eventos que tienen dos posibles resultados, cara-ceca, si-no, apagado-prendido, sano-enfermo, etc. La distribución binomial la vimos dos capítulos atrás y tiene la siguiente forma:
$$p(x \mid n,p) = {\binom {n}{k}}p^x(1-p)^{n-x}$$Hay que tener cuidado que $p$ en la distribución binomial NO es el \text{valor p} del que veníamos hablando! Si no un parámetro de esta distribución.
Según el paradigma frecuentista el procedimiento sería:
In [4]:
p = 0.5
dist_monedas = stats.binom(n, p);
x = np.arange(0, n + 1)
y = dist_monedas.pmf(x)
valor_p = dist_monedas.cdf(90) * 2 # por simetría
T = x / n
T_obs = z / n
valor_p = np.sum(y[T <= T_obs]) * 2
plt.vlines(T, 0, y, 'C0')
plt.vlines(T[T <= T_obs], 0, y[T <= T_obs], 'C1',
label='p = {:.3f}'.format(valor_p))
plt.vlines(T[T >= 1-T_obs], 0, y[T >= 1-T_obs], 'C1')
plt.xlim(0.35, 0.65)
plt.yticks([])
plt.legend()
plt.xlabel(r'$T_{\mathcal{\tilde D}}$', fontsize=14);
Esta aproximación tiene varios problemas
Los $\text{valores p}$ NO son:
Gran parte de los problemas derivados de los $\text{valores p}$, provienen de usarlos como parte de la maquinaria de "prueba de hipótesis nula". La otra parte del problema de usar los $\text{valores p}$ proviene quizá de que estos no derivan de un sistema formal de cálculo de probabilidades si no que fueron introducidos de forma totalmente ad hoc. Como veremos más adelante la Estadística Bayesiana provee de algo llamado factores de Bayes que curiosamente es todo lo que usted siempre quiso que un $\text{valor p}$ fuera, pero que el $\text{valor p}$ NO puede ser (aunque tampoco está excento de peros).
Una forma de usar los $\text{valores p}$ sin meter la pata es no pedirle que ofrezcan la información que no pueden ofrecer. Un $\text{valor p}$ debería ser usado con el fin que originalmente se postuló: una forma aproximada, cuando se cuenta con poca información sobre un problema, para intentar estimar (de forma objetiva) si vale la pena seguir mirando los datos o haciendo experimentos a fin de replicar o extender resultados. Entonces en general valores "bajos" de los $\text{valores p}$ indicarían que es posible que tengamos algo interesante. Por ejemplo algunos autores recomiendan de forma muuuy general interpretar los $\text{valores p}$ de la siguiente manera
un p-valor es una forma de cuantificar cuan sorprendente son los datos si asumimos que no hay efecto. No permiten hacer enunciados sobre hipótesis o teorías.
En estadística se le llama verosimilitud (o likelihood) a la expresión:
$$p(x \mid \theta) \tag{1}$$Donde $\theta$ es un parámetro de un modelo, el cual debemos averiguar, y $x$ son los datos. Siguiendo el ejemplo de la moneda $\theta$ podría ser la cantidad de caras sobre el número de tiradas, es decir la fracción de cara obtenidas, $\theta = \frac{z}{n}$.
Una aclaración: desde la perspectiva frecuentista el likelihood no es una densidad de probabilidad (aunque es proporcional a una), por lo que es común en muchos textos el uso de la siguiente notación para hablar del likelihood:
$$\mathcal{L}(\theta; x)$$A los fines prácticos vamos a obviar esta distinción. Lo que nos interesa es averiguar el valor de $\theta$, es importante notar que para cada valor de $\theta$ que usemos en 1 obtendremos un resultado distinto, es decir existen infinitos valores de $\theta$ de donde elegir ¿Cual elegir? Una respuesta razonable es aquel que maximice la expresión 1, es decir el valor de $\theta$ que me da la máxima probabilidad de observar los datos. Este procedimiento se conoce como estimación por máxima verosimilitud o en inglés maximum likelihood estimation (MLE).
Ahora bien, seguimos con un problema la expresión 1 es una expresión general, para poder calcular algo tenemos que elegir una distribución de probabilidad. Como ya vimos al calcular los $\text{valores p}$ podemos usar una distribución binomial para representar la tirada de monedas. Cuando $n=1$ la distribución binomial se convierte en la distribución de Bernoulli, la cual podemos escribir como:
$$p(x \mid \theta) = \theta^x(1-\theta)^{1-x}\tag{2}$$donde como ya vimos $x \in \{0, 1\}$ (ceca o cara). Esto nos estaría diciendo cual es la probabilidad de ver una cara o una ceca dado un valor de $\theta$ al arrojar una moneda (una sola vez). Como dijimos $\theta = \frac{z}{n}$ es un número entre 0 y 1 que representa la cantidad de caras sobre la cantidad total de tiradas (caras + cecas).
Como no tenemos una moneda si no varias quizá queda más claro escribir
$$p(x_1, x_2 ... x_n \mid \theta) = \theta^{x_1}(1-\theta)^{1-x_1} \theta^{x_2}(1-\theta)^{1-x_2} \ldots \theta^{x_n}(1-\theta)^{1-x_n}$$Fijense que como cada tirada es independiente de la anterior el término a la derecha es simplemente una sucesión de productos
Por lo que de forma más sintética podemos escribir:
$$p(x_1, x_2 ... x_n \mid \theta) = \prod_{i=1}^n \theta^{x_i}(1-\theta)^{1-x_i}$$En la práctica se suele usar el logaritmo del likelihood (log-likelihood) y no el likelihood, ya que estos previene el underflow que podría ocurrir al multiplicar sucesivas veces probabilidades (números entre 0 y 1). Usar el log-likelihood suele también simplificar los cálculos (multiplicaciones se vuelven sumas, potencias productos y divisiones restas). Este truco computacional no cambia los resultados ya que el valor de $\theta$ que maximiza al likelihood es el mismo que el que maxima al log-likelihood.
$$\log p(x_1, x_2 ... x_n \mid \theta) = \sum_{i=1}^n log(\theta^{x_i}(1-\theta)^{1-x_i})$$Resumiedo:
El único ingrediente que nos está faltando es una implementación concreta que nos permita obtener números concretos. Hay varias formas de hacer esto, se puede resolver de forma analítica o numérica. Empecemos por uno de los varios métodos numéricos disponibles.
Un método numérico de fuerza bruta y muy simple es el conocido como método de la grilla. Este método consiste en evaluar una función (en este caso el likelihood) en varios puntos de una grilla (en general, pero no necesariamente, equidistantes).
En la siguiente celda está la solución en Python. Las tres primeras lineas son la clave, mientras el resto sirven para graficar la solución. La primer linea genera la grilla de valores entre 0 y 1. La segunda linea devuelve una lista donde cada elemento es la suma de los logpmf (para cada cara o ceca en monedas_obs), para un valor de $\theta$ (definido en la primer linea). La tercer linea devuelve el valor de $\theta$ para el cual llk es máximo.
In [5]:
θ = np.linspace(0, 1, 100)
llk = [stats.bernoulli(p=θ_i).logpmf(monedas_obs).sum() for θ_i in θ]
θ_hat = θ[np.argmax(llk)]
plt.plot(θ, llk)
plt.plot(θ[np.argmax(llk)], np.max(llk), 'bo',
label='$\hat\\theta ={:.2f}$'.format(θ_hat))
plt.yticks([])
plt.ylabel('llk', fontsize=16, rotation=0, labelpad=15)
plt.xlabel('$\\theta$', fontsize=16)
plt.legend();
In [6]:
θ = np.random.uniform(0, 1, 100)
θ.sort()
llk = [stats.bernoulli(p=θ_i).logpmf(monedas_obs).sum() for θ_i in θ]
θ_hat = θ[np.argmax(llk)]
plt.plot(θ, llk)
plt.plot(θ_hat, np.max(llk), 'bo',
label='$\hat\\theta ={:.2f}$'.format(θ_hat))
plt.yticks([])
plt.ylabel('llk', fontsize=16, rotation=0, labelpad=15)
plt.xlabel('$\\theta$', fontsize=16)
plt.legend();
Una versión más avanza de los métodos de Monte Carlo son la familia de métodos conocidos como Markov Chain Monte Carlo (MCMC). Estos últimos son muy usados en muchas disciplinas científicas e ingenieriles y tienen un rol central en estadística Bayesiana.
La razón principal de la utilidad de los MCMC es que el método de la grilla (al igual que el método de MonteCarlo a secas) falla a medida que la cantidad de variables aumenta.
Quienes tengan interés en indagar más sobre este tema pueden jugar con este notebook.
In [7]:
from scipy.optimize import minimize
Como su nombre lo indica minimize minimiza funciones. Las funciones a minimizar son funciones escritas en Python. Recuerden que maximizar una función $f(\cdot)$ equivale a minimizar $-f(\cdot)$.
In [8]:
def neg_llk(θ):
return -np.sum(stats.bernoulli(p=θ).logpmf(monedas_obs))
θ_ini = 0.5
minimize(neg_llk, θ_ini, method='Nelder-Mead')['x']
Out[8]:
Encontrar el valor de $\theta$ que maximiza el likelihood equivale a derivar el likelihood respecto de $\theta$ e igualar a 0. Al hacer esto se encuentra que:
$$\hat \theta = \frac{\sum_{i=1}^n x_i}{n}$$Es decir el número de caras (éxitos) dividido el número total de tiradas (experimentos). Por lo tanto para nuestro problema $\hat \theta$ es:
In [9]:
z / n
Out[9]:
o lo que es equivalente:
In [10]:
np.mean(monedas_obs)
Out[10]:
Resultado que, para nuestra tranquilidad, concuerda con los valores computados con todos los métodos anteriores.
$\hat \theta$ es un ejemplo de estimación puntual, un número que usamos como estimación de un parámetro de una población. En este caso particular la media de la muestra como aproximación de la media de la población.
Suele ser buena idea acompañar la estimación puntual con una medida de la incertidumbre asociada a esa estimación. Una forma de medir esta incerteza es calculando lo que se conoce como intervalo de confianza.
Este intervalo se puede construir a partir de la distribución de muestreo. A medida que $n$ aumenta el teorema del límite central garantiza que dada una población con media finita $\mu$ y varianza finita (y distinta de cero) $\sigma^2$:
$$\hat \mu \sim \mathcal{N} \left(\mu, \frac{\sigma^2}{n}\right)$$Entonces podemos definir un intervalo de confianza en términos de la estimación puntual y el error estándar (SE) de esa media:
$$IC = [\hat \theta - SE, \hat \theta + SE]$$Este sería un intervalo de confianza de $\pm 1 SE$, si quisiéramos ampliar o achicar el rango deberíamos usar $x$ errores estándar, es decir:
$$IC = [\hat \theta - x SE, \hat \theta + x SE]$$En la práctica es común definir los $IC$ en términos de porcentajes y no $SE$, de hecho el porcentaje más comúnmente usado es 95%. Convertir entre porcentajes y $SE$ es simple si asumimos que la distribución es Gaussiana. En ese caso es conocido que un intervalo del 95% estará de forma aproximada dentro de 1.96 errores estándar.
Podemos usar el siguiente código para pasar de porcentajes a cantidades de $SE$
In [11]:
def IC_to_a(ic=95):
a = (100 - ic) / 100
gaussian = stats.norm(0, 1)
l, u = gaussian.ppf([(a / 2), (1 - a / 2)])
x = np.linspace(-4, 4, 100)
y = gaussian.pdf(x)
plt.plot(x, y);
plt.fill_between(x[x < l], y[x < l], alpha=0.5, color='C2')
plt.fill_between(x[x > u], y[x > u], alpha=0.5, color='C2',
label='$\\alpha$ = {:.3f}'.format(a))
plt.vlines([l, u], 0, gaussian.pdf([l,u]),
label='{:.3f} {:.3f}'.format(l, u))
plt.legend(fontsize=14)
plt.draw()
In [12]:
interact(IC_to_a, ic=(0, 100, 5));
Dada la definición de $IC$ y el hecho que usando $\pm 1.96$ errores estándar podemos definir un intervalo de confianza del 95%, usemos Python para calcularlo.
In [13]:
θ_hat = z / n
dist_monedas = stats.bernoulli(p=θ_hat)
std_error = dist_monedas.std() / n**0.5
r = std_error * 1.96
plt.errorbar(x=θ_hat, y=0,
xerr=r,
fmt='o',
label='$\hat \\theta = {:.2f} \pm {:.2f}$'.format(θ_hat, r))
plt.yticks([])
plt.legend();
Como ya dijimos es muy común usar un intervalo con un nivel de confianza del 95%. Pero este valor es arbitrario, tan arbitrario como el valor de $\alpha$ en las pruebas de hipótesis. De hecho ambos valores están relacionados.
Esencialmente un intervalo de confianza contiene a todos los valores de $\hat \theta$ para los cuales $H_0$ NO sería rechazada, para un dado valor de $\alpha$. Los intervalos de confianza se suelen expresar como $(1 - \alpha) 100 %$, por lo que un valor de $\alpha=0.05$ se relaciona con un intervalo de confianza de 95%.
Como sucede con los $\text{valores p}$ resulta que es fácil interpretar los $IC$ de forma errónea. Por ello veamos cual es la interpretación adecuada.
Si repetimos un experimento (infinitas veces) y cada vez calculamos un intervalo de confianza, entonces obtendremos que el X% de esos intervalos contendrán el valor real del parámetro de interés. Es decir el nivel de confianza, digamos el 95%, NO es sobre un intervalo en particular si no que es un intervalo sobre todo el universo posible de intervalos que se podrían calcular con muestras similares a la nuestra. Por lo tanto NO es posible decir que tenemos una confianza del 95% que NUESTRO intervalo contenga el parámetro, en sentido estricto nuestro intervalo o contiene o no contiene al parámetro real. Si quisiéramos hacer enunciados como "tengo un 95% de confianza que el parámetro real se encuentre en tal o cual rango", deberíamos hacer uso de estadística Bayesiana!
La estadística frecuentista tuvo un gran desarrollo a principios del siglo XX. En tiempos donde el poder computacional era fundamentalmente tracción a sangre, es por ello que gran parte de los métodos frecuentistas consisten en una serie de aproximaciones teóricas que permiten simplificar los cálculos.
La amplia disponibilidad de computadoras rápidas y baratas no solo ha simplificado la aplicación de métodos tradicionales, si no que ha permitido el desarrollo de nuevos métodos basados en aproximaciones numéricas que antes eran inaplicables en la práctica o incluso completamente impensadas. Estos métodos computacionales incluyen el remuestreo (frecuentista), y prácticamente toda la estadística Bayesiana moderna y el Machine Learning.
Empecemos por los métodos de remuestreo. Estos métodos se basan en generar submuestras a partir de nuestros datos, cada muestra la usamos para calcular alguna cantidad o ajustar un modelo. El objetivo es obtener información adicional de los datos o modelos.
Dos métodos de remuestreo muy usados son:
El bootstrapping es un método de remuestreo que permite aproximar la distribución de muestreo. Su principal ventaja es su simplicidad lo que permite usarla con estimadores sencillos como los que hemos usado hasta ahora y también con estimadores más complejos.
La idea del bootstrapping consiste en generar un número grande de muestras a partir de (re)muestrear los propios datos. El muestreo se hace con reemplazo:
De esta forma habré logrado calcular miles de veces el estimador de interés por lo que no solo tendré una estimación puntual, si no una distribución de valores para el estimador. Usando Python se puede escribir de forma más breve, incluso,que usando español ;-)
In [14]:
rep = 10000
alpha = 0.05
θ_hats = np.array([np.random.choice(monedas_obs, replace=True, size=len(monedas_obs)).mean()
for _ in range(rep)])
θ_hat = np.mean(θ_hats)
sd = np.std(θ_hats)
θ_hats.sort()
IC = θ_hats[[int(rep * alpha / 2), int(rep*(1 - alpha / 2))]]
plt.yticks([])
plt.hist(θ_hats, bins=20, density=True,
label='$\hat \\theta$ = {:.2f} $\pm$ {:.2f} \n IC={:.2f}-{:.2f}'.format(θ_hat, sd, *IC))
plt.legend();
"Boot-straps" son las tiras que suelen tener las botas y cuyo propósito es facilitar el colocárselas. Bootstrapping hace referencia a una expresión del inglés que quiere decir algo así como levantarse uno mismo tirando de las tiras de las propias botas, en definitiva algo imposible o que no tiene sentido. De todas formas existen justificaciones teóricas para este procedimiento que garantizan, que bajo ciertas condiciones, las aproximaciones son buenas. Además de justificaciones empíricas que han mostrado la utilidad de este método.
Existen múltiples variantes de este idea general, la que hemos presentado se llama bootstap no-paramétrico. Ya que el método NO asume una distribución para lo datos, simplemente se remuestre. Algunas de las variantes de bootstrapping son:
Bootstrapping paramétrico, a partir de un distribución (por ej la distribución binomial) se obtiene un estimador puntual (por ej usando Maximun Likelihood) y se usa esa distribución (parametrizada con ese estimador puntual) para generar muestras aleatorias. Es decir se ajustan los datos a una distribucion de probabiliad y luego se toman muestras de esa distribución.
Bootstrapping suavizado, lo que se hace es agregar un pequeña cantidad de ruido (generalmente Gaussiano). Esto puede ayudar a llenar los espacios entre datos. De esta forma se obtiene una distribución de muestreo que es más suave. En cierta forma equivale a muestrear no de los datos si no de la Kernel density estimation de los datos. Además este método se puede considerar como un intermedio entre el bootstrapping no-paramétrico y el paramétrico.
Bootstrapping de a bloques, se usa con datos correlacionados por ejemplo series temporales en la que un simple bootstrapping destruiría la correlación existente en los datos.
En general se enseña (¿vende?) al bootstrapping como un procedimiento o método que no asume nada sobre los datos. Es posible mostrar que el bootstrapping es en realidad un modelo estadístico que asume ciertas cosas y que es de hecho un caso especial de modelo Bayesiano (esto lo discutiremos en detalle en el tercer curso de esta especialidad). El creador del Bootstrapping se refiere al mismo como un "poor's man Bayesian method", algo así como la versión "trucha/berreta" de un método Bayesiano. Este método en realidad asume que:
Aún cuando la mayoría de los conjuntos de datos satisfacen estos dos criterios, en la práctica el método (modelo) de bootstrapping suele ser una buena aproximación siempre y cuando se cumplan condiciones como:
Es decir criterios que uno espera, que en general, tenga la mayoría de las muestras con las que uno querría trabajar.
La regresión lineal es uno de los modelos más simples y útiles en el análisis de datos. Es útil en si mismo y como base para otros modelos.
Tenemos una regresión lineal cuando queremos modelar la relación entre una variable $x$ y una variable $y$ usando una linea recta. En el caso más simple, llamado regresión lineal simple, tanto $x$ como $y$ son vectores representando pares de variables. altura-peso, temperatura-humedad, número de premios nobel por país-consumo de chocolate per capita, etc. Una forma de representar este modelo es
$$y_i = \alpha + \beta x_i $$En este modelo $x$ se conoce como variable independiente (o variable predictora) e $y$ se conoce como variable dependiente (o variable predicha). $\beta$ es la pendiente e indica el cambio en $y$ por unidad de cambio de $x$, $\alpha$ indica el valor de $y_i$ para cuando $x_i=0$ y gráficamente representa el punto en que la linea intercepta a la ordenada al origen.
Una forma de resolver este problema es plantearlo de forma probabilista:
$$\mathbf{y} \sim \mathcal{N}(\mu=\alpha + \mathbf{x}\beta, \sigma=\epsilon) $$Es decir $y$ se distribuye como una Gaussiana centrada en $\alpha + \mathbf{x}\beta$, con desviación estándar $\epsilon$.
Veamos como resolver este problema usando el método de máxima verosimilitud y bootstrapping. Para ello generemos unos datos sintéticos.
In [15]:
np.random.seed(42)
x = np.random.normal(0, 1.5, 100)
y = np.random.normal(x, 1)
plt.plot(x, y, 'o');
El likelihood en este caso viene dado por una distribución Gaussiana y ahora nuestras incógnitas son $\alpha$, $\beta$ y $\sigma$. En el siguiente código representadas por una lista/array de tres elementos. El resto es esencialmente lo mismo que hicimos con la moneda.
In [16]:
def neg_llk_r(θr, x, y):
mu = θr[0] + θr[1] * x
return -np.sum(stats.norm(mu, θr[2]).logpdf(y))
θr_ini = [1., 1., 1.]
pred = minimize(neg_llk_r, θr_ini, args=(x, y), method='Nelder-Mead')['x']
plt.plot(x, y, 'o')
plt.plot(x, pred[0] + pred[1] * x, 'k');
Además de encontrar la mejor linea recta que ajusta los puntos podríamos querer tener una medida de la incertidumbre que tenemos de esa recta, eso es fácil de calcular usando usando bootstrapping.
In [17]:
def neg_llk_r(θr, x, y):
mu = θr[0] + θr[1] * x
return -np.sum(stats.norm(mu, θr[2]).logpdf(y))
XY = np.stack((x, y), 1)
θr_ini = [1., 1., 1.]
len_XY = len(XY)
for i in range(0, 50):
XY_bs = XY[np.random.randint(0, len_XY, len_XY)]
x_ = XY_bs[:,0]
y_ = XY_bs[:,1]
pred_bs = minimize(neg_llk_r, θr_ini, args=(x_, y_), method='Nelder-Mead')['x']
plt.plot(x, pred_bs[0] + pred_bs[1] * x, '0.7')
plt.plot(x, y, 'o')
plt.plot(x, pred[0] + pred[1] * x, 'k');
El método es un poco lento, por que tenemos que minimizar cada vez que remuestreamos. El código de la celda anterior es tan solo un ejemplo didáctico y no pretende ser eficiente. Pueden probar stats.linregress en vez del método de optimización.
Así como es posible calcular y visualizar la incerteza para la recta que ajusta los puntos es posible hacer lo mismo para los valores de la variable $y$.
In [18]:
def neg_llkr(θr):
mu = θr[0] + θr[1] * x
return -np.sum(stats.norm(mu, θr[2]).logpdf(y))
θr_ini = [1., 1., 1.]
pred = minimize(neg_llkr, θr_ini, method='Nelder-Mead')['x']
y_pred = pred[0] + pred[1] * x
for i in range(0, 100):
y_bs = np.random.normal(y_pred, pred[2])
plt.plot(x, y_bs, 'ko', alpha=0.02)
plt.plot(x, pred[0] + pred[1] * x, 'C1');
Calcular el $\text{valor-p}$ asumiendo como hipotesis nula que no hay diferencia entre las medias de sepal_width para virginica y versicolor. Asumir que los datos se distribuyen como Gaussianas. Tomar la media y desviación estándar de los datos.
Que sucede con el resultado del punto anterior si la cantidad de datos, $n$, fuese 10? y si fuese 70? Antes de modificar el código y repetir la simulación piense que espera ver. Una vez repetidas las simulaciones para los distintos $n$ Los resultados son lo que esperaba? Qué opina del resultado?
Estime el valor de la media y desviación estándard del sepal_lenght de virginica usando máxima verosimilitud (Método de la grilla y Optimización numérica). Compare ambos valores.
Utilizando la ley de los grandes números estime el intervalo de confianza del 95% para la media de del sepal_lenght de virginica (utilize el valor de $\hat \theta$ estimado por optimización numérica).
Usando bootstrapping estime la media del sepal_lenght de virginica y el intervalo de confianza de 95%. Compare con los valores obtenidos anteriormente.
Modifique el ejemplo anterior para obtener:
Usando el código visto en este notebook para la regresión lineal simple. Comparar las variables petal_lenght vs petal_width para:
Graficar la recta que mejor ajusta y su incertidumbre asociada. En la leyenda incluya los valores coefficiente de correlación de Pearson y la ecuación de la recta (con los valores estimados de $\alpha$ y $\beta$.