NumPy (Numerical Python) es un módulo de Python para cómputo científico. Esta biblioteca contiene muchas funciones útiles en computación científica entre las que se encuentran manipulación de arreglos numéricos, operaciones de álgebra lineal y generación de números pseudo-aleatorios.
Para poder usar NumPy debemos importarlo, la forma más común de importar NumPy es asignándole el alias np:
In [37]:
import numpy as np
NumPy usa una estructura de datos llamada arreglos (arrays). Los arreglos de NumPy son similares a las listas de Python, pero son más eficientes para realizar tareas numéricas. La eficiencia deriva de las siguientes características:
Las listas de Python son muy generales, pudiendo contener objetos de distinto tipo. Los arreglos de NumPy son homogéneos solo pueden contener objetos de un mismo tipo.
En una lista de Python los objetos son asignados dinamicamente, es decir el tamaño de una lista no está predefinidos, siempre podemos agregar más y más elementos. Por el contrario los arreglos de NumPy son estáticos.
Estos dos primeros puntos permiten hacer uso eficiente de la memoria
Otra ventaja de los arreglos es que se comportan de forma similar a los vectores y matrices usados en matemática. Esto facilita muchas de las tareas científicas, precisamente por que el álgebra lineal es el lenguaje usado para pensar y resolver muchos problemas científicos de forma eficiente.
Existen varias rutinas para crear arreglos de NumPy a partir de:
Para crear arreglos a partir de listas (o tuplas) podemos usar la funcion array:
In [38]:
v = np.array([1, 2, 3, 4 , 5, 6])
v
Out[38]:
In [39]:
M = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
M
Out[39]:
El primer arreglo, v, lo creamos a partir de una lista y es por lo tanto unidimensional, mientras que el segundo M lo creamos a partir de una lista anidada (una lista de listas) y resulta en un arreglo bidimensional.
Los arreglos tienen atributos como por ejemplo shape:
In [40]:
v.shape, M.shape
Out[40]:
El shape nos indica la candidad de elemenos en cada eje (o axis). En 2 dimensiones podemos pensarlo como el número de (filas, columnas).
Tambien podemos preguntarle a un array cual es su dimensión:
In [41]:
v.ndim # v es unidimensional
Out[41]:
In [42]:
M.ndim # M es bidimensional
Out[42]:
In [43]:
np.arange(0, 10, 1) # desde, hasta(sin incluir), paso (el paso es opcional!)
Out[43]:
Otra función para crear rangos es linspace que devuelve numeros igualmente espacios en el intervalo [desde, hasta] (es decir incluyendo el hasta). Otra diferencia con arange es que no se especifica el paso sino la cantidad total de números que contendrá el arreglo.
In [44]:
np.linspace(1, 10, 25) # desde, hasta, elementos (elementos es opcional)
Out[44]:
Los números aleatorios son usados en muchos problemas científicos. En la práctica las computadoras son solo capaces de generar números pseudo-aleatorios. Python usa un algortimo llamado Mersenne Twister para generar números pseudo-aleatorios. Este algorítmo es más que suficiente para fines científicos, pero no es útil en caso que necesitemos números pseudo-aleatorios para usar en criptografía. Resumiendo para nuestros fines podemos asumir que los números pseudo-aleatorios que generaremos a lo largo del curso son realmente números aleatorios.
Todas las rutinas de NumPy para generar números aleatorios viven dentro del módulo random.
La función mas simple es rand. Esta función crea un arreglo de números en el intervalo [0, 1). Dentro de ese intervalo los números son equiprobables, es decir es una distribución uniforme. El argumento de rand son las dimensiones del arreglo resultante.
In [45]:
np.random.rand(2, 5) # arreglo con forma (2, 5)
Out[45]:
De forma similar randn (noten la n al final) devuelve muestras a partir de la distribución normal estándar (media = 0, desviación estándar = 1), según las dimensiones que especifiquemos.
In [46]:
np.random.randn(10)
Out[46]:
En Python es común crear una lista vacía que luego se llena de elementos en un loop. En NumPy es común crear un arreglo del tamaño necesario y luego llenarlo de valores. Para estas situaciones resulta conveniente tener a mano funciones que crean arreglos con ceros o unos.
In [47]:
np.zeros((2, 4))
Out[47]:
In [48]:
np.ones((2, 4))
Out[48]:
In [49]:
np.zeros_like(M) # noten que los ceros no tienen ".", es decir son enteros, ¿Por qué?
Out[49]:
In [50]:
np.full((2, 4), 42.)
Out[50]:
Dado que en Python 0 evalua como False y 1 como True, una forma de crear un arreglo de booleanos es:
In [51]:
np.ones((2, 4), dtype=bool)
Out[51]:
In [52]:
datos = np.genfromtxt('datos/microbioma.csv', delimiter=','
, skip_header=1, usecols=(1,2,3,4), dtype='int')
In [53]:
datos.shape
Out[53]:
In [54]:
datos[:4]
Out[54]:
Leamos el archivo microbioma.csv pero esta vez pasando menos argumentos que en el caso anterior:
np.genfromtxt('datos/microbioma.csv', delimiter=',')
¿Cuál son las diferencias entre ambos arreglos?
¿Cómo se explican los nan?
Tip: al menos bajo Linux, es posible usar las celdas de una notebook para ejecutar comandos como si fuese una terminal (o linea de comandos), por ejemplo podemos ejecutar el comando head para ver en encabezado de un archivo.
In [55]:
!head -4 datos/microbioma.csv
In [56]:
M[0] # el primer elemento de M
Out[56]:
La sintaxis usada para indexar y rebanar es una generalización de la usada para las listas de Python. Esta generalización facilita trabajar con arreglos de más de 1 dimensión. Trabajar con más de una o dos dimensiones puede ser un poco confuso, sobre todo al principio ¡Aunque no necesariamente solo al principio!
La siguiente expresión es válida para arreglos
In [57]:
M[0, 1] # el elemento (0, 1) de M
# o también la primer fila de M y de ella el segundo elemento
Out[57]:
Sin embargo esta expresión no es válida para listas.
¿Cuál es la expresión equivalente que funciona con listas?
En la siguiente celda tenemos un ejemplo de una expresión que es común a listas y arreglos.
In [58]:
M[1:] # a partir de la fila 1 todo
Out[58]:
Y este es un ejemplo de rebanado que funciona con arreglos, pero no con listas.
In [59]:
M[1,:] # solo la fila 1 (o la fila 1 en el primer axis y todo en el resto de los axis)
Out[59]:
El poder rebanar/indexar en varias dimensiones en simultaneo nos da flexibilidad para trabajar con subconjuntos de datos contenidos en un arreglo. Veamos más ejemplos:
In [60]:
M[:,1] # solo la columna 1 (o todo en la primer dimensión y la columna 1 en el resto de las dimensiones)
Out[60]:
In [61]:
M[:,1:] # todas las filas y todas las columnas a partir de la columna 1
Out[61]:
In [62]:
M[::-1] # los elementos de M en reversa
Out[62]:
Es importante notar que al tomar rebanadas NumPy NO genera un nuevo arreglo, sino una vista (view) del arreglo original. Por lo tanto si a una rebanada le asignamos un número, se lo estaremos asignando al arreglo original, como se puede ver en el siguiente ejemplo.
In [63]:
M[0, 0] = 42
M
Out[63]:
Distinto es asignar la rebanada a una variable y luego modificar esa variable:
In [64]:
a = M[0, 0]
a = 1
M
Out[64]:
Para crear copias se puede usar la función np.copy() o el método .copy().
Genera una copia de M, llamada K, modifica K y comprueba que M no se modificó
In [65]:
res = np.zeros_like(M, dtype=float)
for i in range(M.shape[0]):
for j in range(M.shape[1]):
res[i, j] = M[i, j] ** 0.5
res
Out[65]:
Otra opción sería usar enumerate:
In [66]:
res = np.zeros_like(M, dtype=float)
for i, fila_i in enumerate(M):
for j, elemento_ij in enumerate(fila_i):
res[i][j] = elemento_ij ** 0.5
res
Out[66]:
NumPy permite vectorizar estas operaciones, es decir podemos calcular la raíz cuadrada de todos los elementos de un arreglo en una sola operación que se aplica a cada uno de los elementos del arreglo:
In [67]:
M**0.5
Out[67]:
Como se ve en el ejemplo anterior vectorizar permite omitir. Esta capacidad de
vectorizar código no está restringida a operadores matemáticos como la potenciación **, funciona con otros operadores y con una gran cantidad de funciones. Por ejemplo
In [68]:
np.sqrt(M)
Out[68]:
In [69]:
np.log(M) # 0j0 log es logaritmo natural
Out[69]:
Funciones como sqrt o log, que operan sobre arreglos elemento-a-elemento se conocen como funciones universales (usualmente abreviadas como ufunc).
Una de las ventajas de usar ufuncs es que permite escribir código más breve. Otra ventaja es que los cómputos son más rápidos que usando loops de Python. Detrás de escena NumPy realiza las operaciones en un lenguaje como C o Fortran, por lo que hay una ganancia considerable en velocidad, respecto de código en Python puro. Además, el código usado por NumPy es código que suele estar optimizado gracias a los años de labor de programadores y científicos que colaboran con proyectos científicos.
Veamos otro ejemplo, como sumar todos los elementos de un arreglo.
In [70]:
np.sum(M)
Out[70]:
En el ejemplo anterior la suma se hizo sobre todos los números contenidos en el arreglo, sin respetar sus dimensiones. En muchas ocaciones es preferible hacer operaciones sobre alguna dimensión en particular, por ejemplo sumar a lo largo de las columnas:
In [71]:
np.sum(M, axis=0)
Out[71]:
O sumar a lo largo de las filas:
In [72]:
np.sum(M, axis=1)
Out[72]:
Un arreglo tendrá tantos axis como dimensiones.
Un característica que facilita vectorizar código es la capacidad de operar sobre arreglos que no tienen las mismas dimensiones. Esto se llama broadcasting y no es más que un conjunto de reglas que permiten aplicar operaciones binarias (suma, multiplicación etc) a arreglos de distinto tamaño.
Consideremos el siguiente ejemplo donde sumamos dos arreglos, elemento a elemento.
In [73]:
a = np.array([0, 1, 2])
b = np.array([2, 2, 2])
a + b
Out[73]:
Podemos ver que el arreglo b contiene 3 veces el número 2. El broadcasting hace que la siguiente operación también sea válida y brinde el mismo resultado que la celda anterior.
In [74]:
a + np.array(2)
Out[74]:
Incluso la siguiente operación es válida:
In [75]:
M + b
Out[75]:
En ambos casos lo que está sucediendo es como si antes de realizar la suma extendieramos una de las partes hasta que las dimensiones coincidan, por ejemplo repetir 3 veces el número 2 o tres veces el vector b.
El broadcasting no funciona para cualquier par de arreglos. La siguiente operación funciona:
In [76]:
M[1:,:] + b
Out[76]:
mientras que la siguiente dará un error
In [77]:
M + b[:2]
El mensaje de error nos dice que NumPy no sabe como hacer para encajar las dimensiones de estos dos arreglos. Considero que este es un error múy razonable ya que no es del todo claro como NumPy podría hacer la operación que le pedidos que haga, además creo que el error es bastante transparente ¿Opinan igual?
Más detalles sobre broadcasting aqui.
In [78]:
M > 6
Out[78]:
Es muy común usar expresiones como la anterior para obtener, de un array dado, solo el subconjunto de valores que cumplen con cierto criterio, como:
In [79]:
M[M > 6]
Out[79]:
o incluso combinando arreglos, como:
In [80]:
M[a == 2]
Out[80]:
In [81]:
np.mean(v)
Out[81]:
Una forma alternativa es usar el método .mean()
In [82]:
v.mean()
Out[82]:
In [83]:
print('varianza {:.2f}'.format(np.var(v)))
print('desviación estándar {:.2f}'.format(np.std(v)))
Además de la varianza existen otras formas de describir la dispersión de los datos. Una de ellas es el rango. Es decir la diferencia entre el valor más grande y el más chico en un conjunto de datos. Un problema con el rango es que es muy sensible a los valores extremos, después de todo se define como la resta de los dos valores más extremos. Una alternativa es calcular un rango pero truncado, es decir dejando de lado valores hacia ambos extremos. Esto se puede hacer con los cuantiles.
Los cuantiles son puntos de corte que dividen al conjunto de datos en grupos de igual tamaño. Existen varios nombres para los cuantiles según la cantidad de divisiones que nos interesen.
En Python el cálculo de estos estadísticos puede realizarse fácilmente usando funciones predefinidas en NumPy.
In [84]:
x = np.random.normal(0, 1, 100)
'percentiles 25={:.2f}; 50={:.2f}; 75={:.2f}'.format(*(np.percentile(x , [25, 50, 75])))
Out[84]:
Un rango que se calcula usando cuantiles y que es muy usado es el rango intercuartil, el cual se calcula como:
$$IQR = q3 − q1 = p75-p25$$y usando NumPy
In [85]:
np.diff(np.percentile(x , [25, 75]))
Out[85]:
Los gráficos ocupan un lugar central en la estadística moderna y en la ciencia de datos, ya sea en el análisis exploratorio de datos o en procesos de inferencia.
Existen varias bibliotecas para hacer gráficos en Python, Matplolib es una de las más usadas. La forma más común de importarla es:
In [86]:
%matplotlib inline
import matplotlib.pyplot as plt
La primer línea es para decirle a la Notebook que los gráficos queden embebidos en la notebook (si no estuvieramos usando la notebook no escribiríamos esta línea).
La segunda línea es la forma estándar de importar matplotlib.
Veamos como hacer un gráfico sencillo.
In [87]:
x = range(20)
y = [i ** 0.5 for i in x]
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel(r'$\sqrt{x}$', rotation=0, labelpad=15);
plot (ya veremos que los hay de otros tipos), donde graficaremos x vs y.Veamos otro ejemplo:
In [88]:
x = range(20)
y = [i ** 2 for i in x]
z = [i ** 1.8 for i in x]
plt.plot(x, y, label=r'$x^2$')
plt.plot(x, z, label=r'$x^{1.8}$')
plt.xlabel('x', fontsize=16)
plt.ylabel(r'$f(x)$', fontsize=16)
plt.legend(fontsize=16);
Existen muchos tipos de gráficos para representar datos. A continuación veremos cinco representaciones comunes para datos unidimensionales:
En un histograma se representa la frecuencia con la que aparecen los distintos valores en un conjunto de datos. Se utilizan barras contiguas para representar los datos. La superficie (y no la altura) de las barras es proporcional a la frecuencia de datos observados. Los datos son agrupados en bins, y suelen graficarse sin normalizar o normalizados. Normalizar implica que el área total del histograma suma 1. No hay que confundir los histogramas con los gráficos de barras que se utilizan para comparar valores discretos entre grupos, mientras que los histogramas se usan para representar distribuciones continuas.
Los histogramas son sensibles a la cantidad de bins que se usan. Si usamos unos pocos bins no lograremos capturar la estructura de los datos, si usamos demasiados bins no solo estaremos representando la estructura de los datos sino también el ruido. Esto se ve más claramente si probamos valores extremos, en el extremo inferior tedríamos a todos los datos representados con una sola barra, en el superior una barra por cada dato.
In [89]:
np.random.seed(440)
x = np.random.gamma(2, 10, size=1000)
plt.hist(x, bins=50, density=True, cumulative=False); # probá cambiar los bins, y los demás argumentos.
Aprovechando lo que hemos aprendido hasta el momento generemos un gráfico que muestre la diferencia entre media y mediana
In [90]:
plt.hist(x, bins=20)
media = np.mean(x)
mediana = np.median(x)
plt.axvline(media, ymax=.9, c='C1', lw='3', label='{:.2f} media'.format(media))
plt.axvline(mediana, ymax=.9, c='C3', lw='3', label='{:.2f} mediana'.format(mediana))
plt.legend(fontsize=14);
Este tipo de gráfico es similar a un histograma, pero en vez de usar un número discreto de bins para representar una distribución se usa una curva suave y continua.
Un gráfico KDE se dibuja de la siguiente forma: se reemplaza cada dato por una distribución Gaussiana y luego se suman todas las Gaussianas. En vez de una distribución Gaussiana es posible usar otras distribuciones. El nombre genérico para esas distribuciones cuya suma se usa como aproximación de una función es el de kernel. La Gaussiana es uno de los kernels más usado.
De forma análoga a lo que sucede con los bins los KDE son sensibles a un parámetro llamado bandwith. Existen varias heurísticas (reglas empíricas que suelen funcionar bien en la práctica) para ajustar el bandwith de forma automática de acuerdo a los datos.
Es posible usar matplotlib para graficar un kde, pero no existe una función que lo haga de forma automática. Es decir es posible pero requiere de cierto trabajo. Lo mismo sucede con otros tipos de gráficos usados para analizar datos, es por ello que existe una biblioteca llamada Seaborn, la cual no es más que una colección de funciones escritas usando matplotlib.
In [91]:
import seaborn as sns
Usando Seaborn, podemos hacer un kde de forma muy simple
In [92]:
sns.kdeplot(x); #también ver la función sns.distplot()
Uno de los problemas de los KDE, o mejor dicho de muchas de sus implementaciones, es que los bordes o límites de las distribuciones no son tenidos en cuenta. Si observamos con cuidado la figura anterior veremos que la curva azul cubre un rango que va aproximadamente desde -10 hasta 100. Pero si miramos los datos veremos que no hay valores negativos ni por encima de $\approx$ 90, es decir el KDE ¡se está inventando datos! Algunas implementaciones, como la utilizada por PyMC3, aplican correcciones que permiten una mejor aproximación en los bordes de la distribuciones. Al usar paquetes como Seaborn que no tienen en cuenta estas consideraciones una forma algo burda de subsanar este problema es truncar el gráfico en los valores mínimos y máximos de los datos. Como Seaborn está escrito usando Matplotlib, podemos hacer esto, y muchas otras modificaciones, usando Matplotlib.
In [93]:
sns.kdeplot(x)
plt.xlim(x.min(), x.max())
plt.xlabel('$x$')
plt.ylabel('$density$');
Los KDE se pueden utilizar para graficar más de una distribución en el mismo gráfico
In [94]:
sns.kdeplot(x, label='$x$')
sns.kdeplot(x**0.9, label='$x^{0.9}$')
sns.kdeplot(10 + x*0.5, label='$10 + 0.5x$')
plt.xlim(x.min(), x.max());
plt.legend();
Un efecto similar puede lograrse con los histrogramas si solo graficamos el contorno de los mismos (histtype='step') o graficando histogramas semitranspartes(alpha menor a 1). Pero en general el efecto no es tan claro como cuando se usan KDEs.
Intenta reproducir el gráfico anterior pero usando histogramas
Este tipo de gráfico sirve para visualizar un conjunto de datos donde una variable es métrica y las demás son categóricas.
Para visualizarlos podemos usar la función stripplot de seaborn (un gráfico similar es el swarmplot). Los stripplot se suelen graficar agregado un poco de ruido (jitter en inglés) a lo largo del eje de las $x$, esto es simplemente un truco para facilitar la visualización de los puntos, que caso contrario caerían todos en una misma línea ya que las variables categóricas no tienen dispersión.
Puede ser útil en si mismo o puede ser usado superpuesto sobre un boxplot o violinplot (ver más adelante).
In [95]:
y0 = np.random.normal(0, 10, size=42)
y1 = np.random.normal(-1, 10, size=50)
sns.stripplot(data=[y0, y1], jitter=True);
Los gráficos de caja están basados en cuartiles. La caja está delimitada (en sus bordes inferior y superior) por el primer y tercer cuartil, mientras que la línea dentro de la caja es el segundo cuartil (la mediana). Los bigotes pueden indicar varias medidas, por eso es siempre importante leer/escribir la leyenda o texto que acompaña a un boxplot, a veces se usa una desviación estándar, otras veces los percentiles 2 y 98, otras veces es una función del rango intercuartil (como en el siguiente gráfico). Los valores por fuera de los bigotes se suelen considerar como datos aberrantes (ver más adelante).
In [96]:
sns.boxplot(data=[y0, y1]);
Los gráficos de violín son una combinación de gráficos de caja con kde.
In [97]:
sns.violinplot(data=[y0, y1]);
Los datos aberrantes (outliers) son valores que están muy alejados de la mayoría de los valores de una distribución. Los datos aberrantes pueden ser errores de medición, errores al procesar los datos o incluso valores correctos pero inusuales (sobre todo cuando la muestra es pequeña). Siempre es buena idea revisar si nuestros datos contienen datos aberrantes y en muchos casos puede llegar a ser conveniente removerlos. Siempre que se remueve un dato aberrante deberá reportarse que fue removido y explicar cual fue el criterio usado para removerlos. Es importante destacar que la decisión de remover datos aberrantes no debe ser tomada a la ligera. Si un supuesto dato aberrante fuese un valor correcto quizá nos estaría indicando que no comprendemos del todo el fenómeno estudiado y al dejarlo de lado podríamos estar perdiéndo información importante!
Existen varios criterios para identificar datos aberrantes. Dos muy usados son:
El primer criterio suele ser usado para distribuciones que se asemejan a Gaussianas, mientras que el segundo es más general ya que el rango intercuartil es una medida más robusta de la dispersión de una distribución que la desviación estándar. El valor de 1.5 es totalmente arbitrario y un valor que se viene usando desde que esta idea fue propuesta. Si nuestros datos son aproximadamente gaussianos, entonces este criterio excluye menos del 1% de los datos.
Según la desigualdad de Chebyshev (y bajo ciertas condiciones), al menos $1 - \frac{1}{k^2}$ de los valores de una distribución están dentro $k$ desviaciones estándar. Es decir casi todos los valores de una distribución de probabilidad están cerca de la media. Por lo tanto el 75% y el 89% de los valores de una distribución se encuentran dentro de 2 y 3 desviaciones estandar, respectivamente. La desigualdad de Chebyshev indica una cota, para varias distribuciones es posible que los valores se encuentren mucho más concentrados alrededor de la media. Por ejemplo esto sucede con las curvas Gaussianas. Para una curva Gaussiana se cumple la regla 68-95-99,7 es decir el 68 por cierto de los datos se encuentra dentro de 1 desviación estándar, el 95 dentro de 2 y el 99.7 dentro de 3.
Los gráficos que hasta ahora hemos visto sirven para visualizar una variable por vez, (aunque sns.kdeplot soporta la visualización de dos variables). Sin embargo, en muchos casos necesitamos entender la relación entre dos variables. Dos variables están inter-relacionadas, si el conocer el valor de una de ellas provee de información sobre el valor de la otra.
Un gráfico de dispersión es un gráfico que utiliza las coordenadas cartesianas para mostrar los valores de dos variables en simultáneo. Estos gráficos son la forma más simple de visualizar la relación entre dos variables.
Supongamos que tenemos dos variables, que creativamente llamaremos $x$, $y$.
In [98]:
x = np.random.normal(size=1000)
y = np.random.normal(loc=x, scale=1)
Usando matplotlib podemos graficar ambas variables usando la función scatter
In [99]:
plt.scatter(x, y, alpha=1)
plt.xlabel('x')
plt.ylabel('y');
Seaborn provee de múltiples opciones para visualizar relaciones entre dos variables, varias de ellas están contenidas en la función joinplot. Esta función además de mostrar el gráfico de dispersión muestra las distribuciones marginales de $x$ e $y$.
In [100]:
sns.jointplot(x, y, kind='scatter', stat_func=None);
Cuando tenemos una gran cantidad de datos los puntos de un scatter plot empiezan a superponerse y puede que ciertos patrones pasen desapercibidos. En estos casos puede ser conveniente agrupar los datos de alguna forma en vez de mirar los datos en crudo.
El siguiente gráfico usa un kernel density estimation tanto para las distribuciones marginales como para la distribución conjunta.
In [101]:
sns.jointplot(x, y, kind='kde', stat_func=None);
Una alternativa a un KDE bidimensional es el hexbin. Este tipo de gráfico es una versión 2D de un histograma. El nombre se debe a que los datos son agrupados en celdas hexagonales. ¿Por qué hexágonos en vez de cuadrados o triángulos? Simplemente porque las celdas hexagonales introducen una menor distorsión en los datos que otras opciones. Esto se debe a las siguientes razones:
Los hexbin son útiles cuando necesitamos visualizar muchos datos. Por muchos me refiero a números por encima de las centenas de miles de datos. Una ventaja de los hexbin sobre los KDE es que el costo computacional es menor.
In [102]:
sns.jointplot(x, y, kind='hex', stat_func=None); # ver también plt.hexbin();
Una alternativa para evitar que algunos puntos opaquen al resto, en un gráfico de dispersión común, es hacer los puntos semitransparentes. En matplotlib la transparencia de los objetos es controlada mediante un parámetro llamado alpha que va entre 0 y 1. Este es un buen momento para volver algunas celdas atrás y ver como este y otros parámetros pueden ser usados para modificar las gráficas realizadas.
Al trabajar con dos variables es común preguntarse por la relación entre ellas. Decimos que dos variables están relacionadas si una provee información sobre la otra. Si en cambio una variable no ofrece información sobre otra decimos que son independientes.
La correlación es una medida de la dependencia de dos variables. Existen varios coeficientes de correlación el más comunmente usado es el coeficiente de correlación de Pearson. Este coeficiente solo sirve para medir relaciones lineales entre variables. El coeficiente de correlación de Pearson es el resultado de dividir la covarianza de las dos variables ($E[(\bar x - x)(\bar y - y)]$) por el producto de sus desviaciones estándar:
$$\rho_{(x,y)}={E[(\bar x - x)(\bar y - y)] \over \sigma_x\sigma_y}$$En palabras (que puede ser más oscuro que en fórmulas), el coeficiente de correlación de Pearson indica como cambia una variable al cambiar la otra respecto de la variación intrínseca de ambas variables.
¿Por qué usar el coeficiente de Pearson y no directamente la covarianza? Fundamentalmente por que la interpretación es más simple. Al dividir por el producto de las varianzas estámos estandarizando la covarianza, obteniendo un coeficiente que no tiene dimensión y que solo puede variar entre -1 y 1 sin importar las unidades de nuestras variables.
La función joinplot, que vimos en el apartado anterior, por defecto nos devuelve el valor del coeficiente de correlación de Person, junto con un valor p cuyo significado estudiaremos en el capítulo 4.
In [103]:
sns.jointplot(x, y, kind='scatter');
Identificar correlaciones puede ser útil para entender como dos variables se relacionan y para predecir una a partir de la otra. Es por ello que muchas veces además de visualizar la relación entre variables se estiman modelos que ajustan a los datos. Como por ejemplo líneas rectas. En el tercer curso de este programa veremos como crear modelos lineales y no-lineales. Por ahora simplemente nos conformaremos con dejar que seaborn ajuste los datos a un recta por nosotros.
In [104]:
sns.jointplot(x, y, kind='reg');
En la siguiente imagen se puede ver varios conjuntos de datos y sus respectivos coeficientes de correlación de Pearson. Es importante notar que el coeficiente de correlación de Pearson refleja la linearidad y la dirección de dicha linearidad (primera fila), pero no la pendiente de dicha relación (fila del medio). Tampoco es capaz de capturar relaciones no-lineales. En la fila del medio la línea con pendiente cero tiene un coeficiente de correlación de Pearson indefinido, ya que la varianza de la variable $y$ es 0.
<img src='imagenes/Correlación.png' alt="correlación", width=600, height=600>
Si existe algún tipo de mecanismo que hace que una variable dependa de otra deberá existir correlación (aunque no necesariamente lineal). Pero lo opuesto no es siempre cierto, dos variables pueden estar correlacionadas sin que exista ningún tipo de mecanismo que las vincule. Dado el gran conjunto de variables que es posible medir no debería ser sorprendente que existan correlaciones espurias. Por ejemplo en la siguiente figura se puede ver que el número de piratas y la media de la temperatura global están inversamente correlacionados.
<img src='imagenes/pirates_temp.png' alt="Pirates_temp", width=600, height=600>
Este gráfico fue creado a propósito para ilustrar, entre otros puntos, que correlación no implica causalidad (nótese además que el orden de los datos en el eje $x$ es erróneo y la escala es al menos problemática). Para más detalles del origen de esta gráfica leer esta entrada de wikipedia
La aparente relación entre las variables temperatura media y cantidad de piratas podría ser explicada de varias formas, quizá es pura casualidad o quizá se podría establecer que los cambios introducidos por la revolución industrial terminaron por un lado aumentando la cantidad de $CO_2$ (y otros gases de invernadero) y por el otro produciendo cambios socio-culturales y tecnológicos que llevaron (luego de una larga cadena de sucesos) a la disminución de los piratas. Pero ¡no es cierto que podamos contrarrestar el calentamiento global simplemente aumentando la cantidad de piratas!
Para poder establecer una relación causal a partir de una correlación hace falta poder establecer y probar la existencia de un mecanismo que vincule ambas variables. Espero que este ejemplo haya servido para ayudarles a entender que correlación no implica causalidad.
NumPy junto con SciPy son el núcleo de todo el ecosistema científico de Python. Al trabajar en computación científica muchas veces necesitaremos acceso a rutinas numéricas para realizar tareas como interpolar, integrar, optimizar, realizar análisis estadísticos, procesar audio, procesar imágenes, etc.
SciPy es una librería de computación científica construida encima de NumPy que nos ofrece muchas de estas funciones. Al igual que como sucede con NumPy, SciPy cuenta con muchas rutinas rápidas, confiables y facilmente disponibles. SciPy es también el nombre de un grupo de conferencias donde participan usuarios y desarrolladores de herramientas de computación científica en Python.
No es muy común importar todas las funciones de SciPy, en general se importan submodulos o incluso funciones particulares. Si necesitaramos funciones estadísticas probablemente las importaríamos de la siguiente forma:
In [105]:
from scipy import stats
Y luego lo usaríamos como:
In [106]:
stats.describe(x)
Out[106]:
O podríamos querer calcular una regresión lineal:
In [107]:
stats.linregress(x, y)
Out[107]:
En general para usar SciPy no hay que aprender mucho más que lo que ya aprendimos de NumPy y luego leer la documentación de cada función específica que necesitemos usar. En los próximos capítulos veremos algunos ejemplos de su uso a medida que vayamos necesitándolo.
Vamos a explorar que otras funciones nos ofrece NumPy, para ello podemos usar los métodos de instrospección que ofrece Jupyter y la documentación de NumPy·
Existen muchas otras operaciones que se pueden hacer sobre arreglos y muchas otras funciones que ofrece NumPy. Exploremos las siguientes:
reshapeconcatenatehstackvstacksplitflattensortargsortloadtxtCrear un array llamado arr con dimensiones 3x3 y que contenga los enteros de 0 hasta 9.
arrCrear un array como arr pero con los números impares reemplazados por -1
Generar datos gaussianos con np.random.randn(size=s) donde s es igual a 10, 100 o 1000 y para cada caso vamos a contar cuantos puntos son outliers de acuerdo a la regla del rango intercuartil (usando el valor de 1.5) y cuantos valores son aberrantes usando 2 y 3 desviaciones estándar. Para asegurarse de tener números confiables repetir el ejercicio varias veces para cada s y reportar el número promedio de datos aberrantes y su desviación estándar.
Comparar las formulas de la varianza y covarianza y explicar por que es correcto decir que la varianza es la covarianza de una variable respecto de ella misma.
Escribir una función que calcule el coeficiente de correlación de Pearson. Evaluar que el resultado es correcto comparándola con el resultado de stats.linregress
Generar 3 o 4 conjuntos de datos utilizando NumPy. Podés intentar combinando métodos random, funciones trigonométricas, logaritmos, etc.
Hacer gráficos unidimensionales de los datos generados en el punto anterior. Hacer al menos un ejemplo usando solo matplotlib, otro usando solo seaborn y otro usando matplotlib junto con seaborn.
Hace al menos un gráfico bidimensional que refleje la relación entre dos variables generadas en el punto 6. Incluir como parte de la leyenda, el valor del coeficiente de correlación de Pearson (usando tanto la función creada en 5 como stats.linregress)