Tutorial sobre estadística en Python: I Distribuciones, Distribuciones, g ráficos, clustering
In [1]: import import import import
pandas as pd numpy as np matplotlib.pyplot as pl plt t scipy.stats as s
np.__version__ In [2]: print np.__version__ import scipy scipy.__version__ print scipy.__version__ matplotlib.__version__ print matplotlib.__version__ 1.8.0 0.9.0 1.3.1
Referencias Este
es
un
tutorial
sobre
análisis
estadístico
en
Python:
tutorial
R.
Olson
(http://www.randalolson.com/2012/08/06/statis (http://www.randalolson.com/2012/08/06/statistical-analys tical-analysis-made-easy-in-python/) is-made-easy-in-python/) Esta
es
la
documentación
de
la
librería
estadística
de
SciPy:
scipy.stats
(http://docs.scipy.org/doc/s (http://docs.scipy.org/doc/scipy/reference/tutorial/s cipy/reference/tutorial/stats.html) tats.html) Otro tutorial también sobre análisis estadístico y Python, en forma de notebook: Statistical Data Analysis in Python (https://github.com/fonnesbeck/statistical-analysis-python-tutorial), con una parte interesante dedicada a modelos Y un curso con el material disponible en forma de notebooks: Coursera data analysis course, done in Python (https://github.com/herrfz/dataanalysis). Aunque he descargado todos los notebooks en la carpeta Coursera_dataanalysis,
lo
mejor
es
utilizar
el
índice
del
fichero
readme
Data
Analysis
(https://github.com/herrfz/dataanalysis/blob/mast (https://github.com/herrfz/dataanalys is/blob/master/REA er/README.md) DME.md) El
material
completo
del
curso
de
Cursera
está
aquí:
Coursera
(https://github.com/jtleek/dataanalysis) , aunque los pdf de las lecciones los he descargado en Dropbox/personal/R proyectos El blog Micropore (http://micropore.wordpress.com/) tiene bastantes cosas de Python y astronomía
Generación de valores aleatorios con Numpy
In [3]: # # # #
gene g enerac ración ión de val valore ores s alea aleator torios ios seg según ún dis distri tribuc bución ión nor normal mal normal( norm al(loc loc=0. =0.0, 0, sca scale= le=1.0 1.0, , size size=No =None) ne) loc es la med loc media, ia, sca scale le la des desvia viació ción n está estánda ndar r size es un un val size valor or ent entero ero
# o una tup tupla la con el sha shape pe des desead eado o alturas = np.random.normal(188, 3, 20) print alturas [
190.09944441 190.09944441 187.08032995 196.58267576 191.81629275
193.31053145 192.10140399 189.08948296 194.01804315 190.53477337 188.38027873 190.74397206 191.10964306 187.70675913 191.99139412 188.14624532 190.57930883 185.77887953 189.25968599 183.90044742 187.989887 ]
In [4]: # Cuan Cuando do des deseam eamos os val valore ores s segú según n la la dist distrib ribuci ución ón nor normal mal est estánd ándar ar # (mea (mean=0 n=0, , stde stdev v =1), =1), pod podemo emos s util utiliza izar r # stan standar dard_n d_norm ormal( al(siz size), e), don donde de siz size e pued puede e ser ser un val valor or ent entero ero o un un shap shap e np.random.standard_normal((4,2)) Out[4]: array([[-0.12656278, [ 0.96472613, [ 0.40339144, [-1.15241441,
1.04534887], -0.95483344], -0.95483344], -2.29152082], -2.29152082], 0.73026361]])
In [22]: # Gene Generac ración ión de val valore ores s alea aleator torios ios seg según ún dis distri tribuc bución ión bin binomi omial al # bino binomia mial(n l(n, , p, p, size size) ) coinFlips = np.random.binomial(1, 0.5, 10) print coinFlips [0 1 1 1 1 1 0 0 1 1]
In [9]: # # # #
Gene G enerac ración ión de val valore ores s alea aleator torios ios seg según ún la dis distri tribuc bución ión uni unifor forme me En un En un inte interva rvalo lo sem semiab iabier ierto to [lo [low, w, hig high) h) uniform unif orm(lo (low=0 w=0.0, .0, hig high=1 h=1.0, .0, siz size=1 e=1) ) size pue size puede de ser un ent entero ero o una una tup tupla la ind indica icando ndo un sha shape pe
np.random.uniform(size=10) Out[9]: Out[ 9]: arra array([ y([ 0.552651 0.5526518 8 , 0.24854146,
0.59462721, 0.594627 21, 0.7054177 ,
0.13537618, 0.135376 18, 0.17873533,
0.91364061, 0.9136406 1, 0.70354042,
0.59570865, 0.5957086 5, 0.30468929])
In [3]: # # # #
gene g enerac ración ión de val valore ores s alea aleator torios ios seg según ún dis distri tribuc bución ión nor normal mal normal( norm al(loc loc=0. =0.0, 0, sca scale= le=1.0 1.0, , size size=No =None) ne) loc es la med loc media, ia, sca scale le la des desvia viació ción n está estánda ndar r size es un un val size valor or ent entero ero
# o una tup tupla la con el sha shape pe des desead eado o alturas = np.random.normal(188, 3, 20) print alturas [
190.09944441 190.09944441 187.08032995 196.58267576 191.81629275
193.31053145 192.10140399 189.08948296 194.01804315 190.53477337 188.38027873 190.74397206 191.10964306 187.70675913 191.99139412 188.14624532 190.57930883 185.77887953 189.25968599 183.90044742 187.989887 ]
In [4]: # Cuan Cuando do des deseam eamos os val valore ores s segú según n la la dist distrib ribuci ución ón nor normal mal est estánd ándar ar # (mea (mean=0 n=0, , stde stdev v =1), =1), pod podemo emos s util utiliza izar r # stan standar dard_n d_norm ormal( al(siz size), e), don donde de siz size e pued puede e ser ser un val valor or ent entero ero o un un shap shap e np.random.standard_normal((4,2)) Out[4]: array([[-0.12656278, [ 0.96472613, [ 0.40339144, [-1.15241441,
1.04534887], -0.95483344], -0.95483344], -2.29152082], -2.29152082], 0.73026361]])
In [22]: # Gene Generac ración ión de val valore ores s alea aleator torios ios seg según ún dis distri tribuc bución ión bin binomi omial al # bino binomia mial(n l(n, , p, p, size size) ) coinFlips = np.random.binomial(1, 0.5, 10) print coinFlips [0 1 1 1 1 1 0 0 1 1]
In [9]: # # # #
Gene G enerac ración ión de val valore ores s alea aleator torios ios seg según ún la dis distri tribuc bución ión uni unifor forme me En un En un inte interva rvalo lo sem semiab iabier ierto to [lo [low, w, hig high) h) uniform unif orm(lo (low=0 w=0.0, .0, hig high=1 h=1.0, .0, siz size=1 e=1) ) size pue size puede de ser un ent entero ero o una una tup tupla la ind indica icando ndo un sha shape pe
np.random.uniform(size=10) Out[9]: Out[ 9]: arra array([ y([ 0.552651 0.5526518 8 , 0.24854146,
0.59462721, 0.594627 21, 0.7054177 ,
0.13537618, 0.135376 18, 0.17873533,
0.91364061, 0.9136406 1, 0.70354042,
0.59570865, 0.5957086 5, 0.30468929])
In [15]: # Gene Generac ración ión de val valore ores s alea aleator torios ios ent entero eros s segú según n la la dist distrib ribuci ución ón uni unifor for me # en en el el inte interva rvalo lo sem semiab iabier ierto to [lo [low, w, hig high) h) # o bie bien, n, si hig high h se se omit omite, e, en [0, low low) ) # size size es un ent entero ero o una una tup tupla la ind indica icando ndo un sha shape pe #randi #ra ndint( nt(low low, , high high=No =None, ne, siz size=N e=None one) ) np.random.randint(0,12,(3,4)) ,4)) print np.random.randint(0,12,(3 # Tamb También ién se pue puede de uti utiliz lizar ar ran random dom_in _integ tegers ers, , simi similar lar a la la ante anterio rior r # pero pero que gen genera era val valore ores s ente enteros ros en el int interv ervalo alo cer cerrad rado o #[low, #[l ow, hig high] h] print np.random.random_integers(1,6, (1,6, size=10) print np.random.random_integers [[ 5 2 0 [ 8 11 11 [11 9 0
8] 7] 9]]
[6 6 3 4 1 4 1 6 6 1]
Obtención de muestras con Numpy In [16]: # # # # # # #
rand r andom. om.cho choice ice(a, (a, siz size=N e=None one, , repl replace ace=Tr =True, ue, p=N p=None one) ) a es un arr array ay 1-D obl obliga igator toriam iament ente, e, o bien bien un num num. . ente entero ro (ver eje (ver ejempl mplo o sigu siguien iente) te) size es un val size valor or ent entero ero o una una tup tupla la ind indica icando ndo un sha shape pe p es un arr array ay de pro probab babili ilidad dades es de igu igual al lon longit gitud ud que a si se si se omit omite e se se supo supone ne dis distri tribuc bución ión uni unifor forme me par para a a replace repl ace=Fa =False lse si se des desea ea una mue muestr stra a sin sin rem rempla plazam zamien iento to
np.random.choice(alturas np.random.choice(alturas, , size=5, replace=False) replace=False) 196.58267576, 190.74397206, 190.74397206, Out[16]: array([ 196.58267576, 189.08948296])
189.25968599, 189.25968599,
191.99139412, 191.99139412,
In [18]: # Otra Otra pos posibi ibilid lidad ad es hac hacer er a igua igual l a un val valor or ent entero ero n # Vien Viene e bien bien par para a por por eje ejempl mplo o inde indexar xar con él # En En este este cas caso o las las mue muestr stras as se tom toman an del arr array ay ara arange nge(n) (n) np.random.choice(10, np.random.choice(10, size = (2,3), replace=False) Out[18]: array([[1, 2, 8], [3, 9, 4]])
Fijar una semilla
In [71]: # Vamos obteniendo valores diferentes en cada llamada a la función np.random.seed(12345) print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) [3 6 6 2 5 2 6 3 6 2] [2 4 2 4 1 3 2 4 6 3] [2 3 4 6 1 6 2 4 5 1]
In [76]: # Y si volvemos a activar la semilla, recomenzamos la misma secuencia np.random.seed(12345) print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) [3 6 6 2 5 2 6 3 6 2] [2 4 2 4 1 3 2 4 6 3] [2 3 4 6 1 6 2 4 5 1]
In [77]: # Sin embargo una nueva celda del notebook # comienza simepre con una semilla aleatoria # a menos que volvamos a fijar una semilla print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) print np.random.random_integers(1,6, size=10) [4 5 5 6 4 4 2 4 6 3] [6 4 1 6 3 5 6 4 4 5] [4 1 4 1 1 4 6 3 5 5]
Estadística descriptiva Media
In [295]: np.random.seed(12345) x =np.random.randint(0,12,(3,4)) print x # Obtener la media de todos los valores del array print "media total: ", x.mean() # Obtener la media de cada columna print "media calculada variando la fila: ", x.mean(axis=0) # Obtener la media de cada fila print "media calculada variando la columna: ", x.mean(axis=1) [[ 2 [ 9 [ 6 media
5 1 4] 5 2 1] 1 11 9]] total: 4.66666666667
media calculada variando la fila: [ 5.66666667 3.66666667 4.66666667] media calculada variando la columna: [ 3. 4.25 6.75]
Varianza y desviación estándar In [291]: # varianza muestral (dividiendo por n-1) print x.var(ddof=1) print x.var(ddof=1, axis=0) print x.var(ddof=1, axis=1) 12.2424242424 [ 12.33333333 [ 3.33333333
5.33333333 12.91666667
30.33333333 16.33333333] 18.91666667]
In [292]: # Desviación típica muestral print x.std(ddof=1) print x.std(ddof=1, axis=0) print x.std(ddof=1, axis=1) 3.49891758154 [ 3.51188458 2.30940108 [ 1.82574186 3.59397644
5.50757055 4.04145188] 4.34932945]
Relación lineal entre dos variables En este apartado utilizaremos el paquete estadístico de scipy, que se importa con: "import scipy.stats as s"
Coeficiente de correlación r de Pearson
4.66666667
In [54]: # Construyamos dos arrays x e y con los que experimentar np.random.seed(12345) x = np.random.uniform(0,10,size=20) noise = np.random.uniform(0, 3, size=20) y = 5 + x + noise plt.scatter(x,y);
In [45]: r, p = s.pearsonr(x,y) print r, p 0.966159857324 4.88820573698e-12
El coeficiente de determinación es el cuadrado del coeficiente de correlación porcentaje de variabilidad de la variable
, y se interpreta como el
que es explicado por el modelo lineal.
In [39]: # El valor de r**2 es el coeficiente de determinación, # que indica la fortaleza de la relación entre ambas variables
print u"El coeficiente de determinación es %.2f" %r**2 # El 93% de la variación de y se puede explicar por la variable x, # el resto será debido a causas desconocidas, variables ocultas # o variabilidad inherente El coeficiente de determinación es 0.93
Regresión lineal Cálculo de los parámetros de la recta de regresión lineal se hacen con la función linregress() de Scipy. http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mstats.linregress.html
In [48]: slope, intercept, r_value, p_value, std_err = s.linregress(x, y)
In [49]: print r_value 0.966159857324
In [53]: plt.scatter(x,y) plt.plot([0, 10], [intercept, slope * 10 + intercept], '-r', lw=2);
Factores / variables categóricas In [228]: # Los factores de R equivalen a las variables categóricas de Pandas ciudades = np.array(['Madrid', 'Cuenca', 'Toledo', 'Madrid', 'Soria', ' Toledo', 'Madrid']) cF = pd.Categorical.from_array(ciudades) cF Out[228]: Categorical: [Madrid, Cuenca, Toledo, Madrid, Soria, Toledo, Madrid] Levels (4): Index(['Cuenca', 'Madrid', 'Soria', 'Toledo'], dtype=object) In [231]: # Los niveles de la variable categórica son los valores diferentes que toma cF.levels Out[231]: Index([u'Cuenca', u'Madrid', u'Soria', u'Toledo'], dtype=object) In [232]: # Los labels son el contenido de la variable, codificada numéricamente # haciendo referencia a los levels cF.labels Out[232]: array([1, 0, 3, 1, 2, 3, 1])
Las distribuciones de probabilidad con Scipy
Cuestiones de notación
Dada una variable aleatoria valor
, el rango de
,
es el conjunto de valores que puede tomar. Cada
es llamado un cuantil, y la probabilidad de que
tome el valor
se designa por
Un "random variate" es un resultado particular de una variable aleatoria.
Vamos a utilizar el paquete stats de Scipy, que normalmente se importa con: "import scipy.stats as s". Todas las distribuciones de probabilidad tienen en este paquete un tratamiento similar, con los mismos métodos en todas ellas. Por ejemplo, para ver los métodos de la distribución normal, hacer: print s.norm.doc
La distribución normal
In [12]: # pdf: Probability Density Function # s.norm.pdf(array x de cuantiles, loc, scale) # loc es la media, scale es la desviación estándar x = np.linspace(-5, 5, num=100) normalDensity = s.norm.pdf(x, loc=0, scale=1) plt.plot(x, normalDensity) plt.xlabel('cuantiles');
In [13]: # cdf: Cumulative Distribution Function # s.norm.pdf(array de cuantiles, loc, scale) # loc es la media, scale la desviación estándar # Probabilidad de |X| < sigma sigma = 4 print "%.3f %%" %(s.norm.cdf(sigma, loc=0, scale=sigma)-s.norm.cdf(-sig ma, loc=0, scale=sigma)) 0.683 %
In [15]: # ppf: Percent Point Function (Inversa de CDF) s.norm.ppf(0.683/2 + 0.5, loc=0, scale=1) Out[15]: 1.0006418287624492 In [16]: # Devuelve la media y varianza de la distribución s.norm.stats(loc=3, scale=4) Out[16]: (array(3.0), array(16.0)) In [17]: # rvs: generación de valores aleatorios (random variates) # bajo la variable aleatoria considerada ( la normal en el ejemplo) s.norm.rvs(loc=5, scale=2, size=10) Out[17]: array([ 2.07049884, 4.37584035,
5.71139329, 4.60381517,
5.58134749, 7.66738297,
6.19376692, 8.44622194,
5.92728227, 1.63065077])
"Freezing" Con el fin de no tener que escribir en cada caso los parámetros loc y scale de una distribución, se puede definir una variable aleatoria con "rv" del tipo y con los parámetros deseados, y despues referirnos a ella con el nombre de la variable que le hayamos dado. Es decir, el objeto que hemos definido tiene los mismos métodos que la distribución original
In [21]: mi_rv = s.norm(loc=5, scale=2) mi_rv.stats() Out[21]: (array(5.0), array(4.0)) In [24]: mi_rv.rvs(size=10) Out[24]: array([ 6.16432526, 6.35598259,
5.48459398, 4.12032957,
3.25193664, 4.99845849,
Distribución binomial
In [34]: # pmf: Probability Mass Function # Sustituye a pdf en las variables discretas # pmf(x, n, pr) x=[0,1,2] s.binom.pmf(x,2, 0.5) Out[34]: array([ 0.25,
0.5 ,
Gráficos exploratorios Representar muchos puntos
0.25])
3.58073744, 6.55752111,
7.89030591, 3.26517431])
In [297]: x = np.random.normal(size=1e4) y = np.random.normal(size=1e4) plt.scatter(x, y); # con plt.plot(x,y, 'o') hubieramos obtenido el mismo resultado Out[297]:
Aquí no se ve nada, una posibilidad es representar una muestra (consideramos los valores de x e y emparejados)
In [298]: r = np.random.choice(10**4, size=1000, replace=False) plt.scatter(x[r], y[r]);
In [309]: # También puede intentarse variar la transperencia # y tamaño de los puntos plt.scatter(x, y, edgecolors='none', alpha=0.025, s=60);
In [310]: # Y otra opción es utilizar un gráfico de R %load_ext rmagic %Rpush x y %R smoothScatter(x, y, nrpoints=0) The rmagic extension is already loaded. To reload it, use: %reload_ext rmagic
In [8]: # Se puede añadir como opción un color map, por ejemplo: # cmap=plt.cm.Greys # cmap=plt.cm.Reds # # cmap=plt.cm.hot plt.hexbin(x, y) # El color en cada hexágono corresponde a la frecuenci a # También se puede utilizar la opción bins=n para norm alizar cb=plt.colorbar() cb.set_label('frecuencias') fig = plt.gcf() # Get current figure fig.set_size_inches(10,8)
Ahora, el mismo ejemplo utilizando figure y axis explícitamente:
In [9]: fig, ax = plt.subplots() fig.set_size_inches(10,8) im = ax.hexbin(x,y) fig.colorbar(im, ax=ax);
A continuación haremos un histograma 2d, de baja resolución para poder analizarlo
In [40]: fig, ax = plt.subplots() fig.set_size_inches(10,8) hist, xedges, yedges, im = ax.hist2d(x,y, bins=(10,10), range=[[-4.,4.], [-4., 4.]]) fig.colorbar(im, ax=ax);
Vamos a ver lo que significa cada una de las variables devueltas por hist2d:
In [26]: # hist es la frecuencia en cada bin, en este caso un grid de 10x10 print hist.shape print hist (10, [[ [ [ [ [ [ [ [ [ [
10) 0. 0. 0. 0. 2. 4. 0. 1. 1. 0.
0. 0. 4. 10. 17. 19. 9. 4. 0. 0.
0. 7. 26. 62. 134. 130. 93. 24. 6. 0.
1. 18. 69. 234. 445. 456. 239. 90. 12. 0.
2. 17. 136. 449. 845. 824. 451. 137. 13. 0.
5. 19. 151. 455. 823. 824. 435. 151. 20. 3.
0. 10. 61. 220. 477. 447. 253. 78. 15. 2.
0. 6. 16. 67. 135. 147. 72. 20. 6. 0.
0. 0. 6. 9. 27. 23. 14. 3. 0. 0.
0.] 0.] 0.] 2.] 2.] 2.] 1.] 0.] 0.] 0.]]
In [31]: # xbins e ybins son los límites de los intervalos semiabiertos # que definen cada bin print xbins print ybins [-4. [-4.
-3.2 -2.4 -1.6 -0.8 -3.2 -2.4 -1.6 -0.8
0. 0.
0.8 0.8
1.6 1.6
2.4 2.4
3.2 3.2
4. ] 4. ]
Ahora vamos a ver como podemos repr esentar un diagrama de contornos sobre la propia figura. Lo que queremos representar como tercera magnitud son las frecuencias en cada celda... Hay un ejemplo en: http://micropore.wordpress.com/2011/10/01/2d-density-plot-or-2d-histogram/
In [263]: fig, ax = plt.subplots() fig.set_size_inches(10, 8) # Rango de valores en cada eje rango = [[-4., 4.],[-4., 4.]] # Generamos el histograma y el gráfico a la vez # Valores de bins elevados dan mucha fragmentación en las curvas hist, xedges, yedges, im = ax.hist2d(x,y, bins=(30,30), range=rango, cmap=plt.cm.Greys) fig.colorbar(im, ax=ax); extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] # Niveles para las curvas de nivel niveles = [30, 50, 70, 90] colores=['green', 'blue', 'yellow', 'red'] cset = ax.contour(hist, niveles, linewidths=2, extent=extent, colors=co lores) ax.clabel(cset, inline=1, fontsize=15, fmt='%d')
Out[263]:
Pero, supongamos que solo queremos los contornos y no queremos representar la figura del histograma. En ese caso el histograma lo generamos con numpy
In [262]: fig, ax = plt.subplots() fig.set_size_inches(8, 8) # Rango de valores en cada eje rango = [[-4., 4.],[-4., 4.]] # Generamos un histograma wD de frecuencias con numpy hist, xedges, yedges= np.histogram2d(x,y, bins=(30,30), range=rango) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] # Niveles para las curvas de nivel de las frecuencias niveles = [30, 50, 70, 90] colores=['green', 'blue', 'yellow', 'red'] cset = ax.contour(hist, niveles, linewidths=2, extent=extent, colors=co lores) ax.clabel(cset, inline=1, fontsize=15, fmt='%d');
Veamos ahora que significa normalizar los histogramas. Empecemos con el de contornos:
In [200]: rango = [[-4., 4.],[-4., 4.]] hist, xedges, yedges= np.histogram2d(x,y, bins=(30,30), range=rango) total_count = sum(hist) print total_count print hist[14:16,14:16] 9998.0 [[ 111. [ 100.
115.] 111.]]
Ahora normalizamos. El array 2D ahora contiene la densidad en cada bin, es decir, el número de observaciones en el bin dividido por el área del bin.
In [229]: rango = [[-4., 4.],[-4., 4.]] hist, xedges, yedges= np.histogram2d(x,y, bins=(30,30), range=rango) histn, xedges, yedges= np.histogram2d(x,y, bins=(30,30), range=rango, normed=True) print "el resultado es:" print histn[14:16,14:16] # obtenido del siguiente modo: hist_suma_1 = hist/sum(hist) print '\n' area = 8. * 8. / (30*30) print hist_suma_1[14:16,14:16]/area el resultado es: [[ 0.15612497 0.1617511 ] [ 0.14065313 0.15612497]] [[ 0.15612497 [ 0.14065313
0.1617511 ] 0.15612497]]
De hecho, si consideramos que la matriz histn es una matriz de densidades, cumplirá: suma(densidad_i x area bin_i) = 64/900 x suma(densidad_i) = 64/900 * sum(histn) = 1 En efecto:
In [232]: 64./900 * sum(histn) Out[232]: 0.99999999999999833
In [236]: # Este es un ejemplo con valores normalizados. Las curvas dicen # Que por ejemplo la densidad es superior a 0.01 dentro de la curva fig, ax = plt.subplots() fig.set_size_inches(8, 8) # Rango de valores en cada eje rango = [[-4., 4.],[-4., 4.]]
histn, xedges, yedges= np.histogram2d(x,y, bins=(30,30), range=rango, normed=True) extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]] # Niveles para las curvas de nivel de las densidades niveles = [0.01,0.05, 0.1, 0.15] colores=['green', 'blue', 'yellow', 'red'] cset = ax.contour(histn, niveles, linewidths=2, extent=extent, colors=c olores) ax.clabel(cset, inline=1, fontsize=15, fmt='%.2f');
Heatmaps
En lugar de emplear la función plt.hist2d, también se puede generar un histograma 2D con numpy.histogram2dy, y representarlo con imshow():
imshow() se emplea para representar imágenes. En nuestro caso la imagen será una matriz bidimensional, generada con histogram2d(), donde cada celda corresponde a un bin, y en ella hay un valor entero (una cuenta de observaciones que caen dentro del bin). De modo que imshow() lo interpreta como una imagen en escala de grises, aunque lo pinte en color, dependiendo del mapa de color que utilicemos. Esto tiene como consecuencia que imshow() suaviza la imagen (los bins/pixels) no se muestran con claridad.
In [270]: # generamos datos de test con la distribución normal estándar x = np.random.randn(8873) y = np.random.randn(8873) heatmap, xedges, yedges = np.histogram2d(x, y, bins=50) # genera un hea tmap 50x50 extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] plt.imshow(heatmap, extent=extent);
En cambio, matshow() representa la matriz sin hacer ningun suavizado. Por ello matshow es la mejor opción cuando se quiere utilizar un mapa de colores para ver como se distribuyen las observaciones.
In [273]: plt.matshow(heatmap, extent=extent);
Clustering
Agrupamiento jerárquico (hierarchical clustering) La idea del clustering o agrupamiento jerárquico es construir un arbol de smilaridades basado en distancias entre cada dos observaciones.
Referencia
a
la
librería
scipy.cluster.hierarchy
(http://docs.scipy.org/doc/scipy/reference/cluster.hierarchy.html#module-scipy.cluster.hierarchy) Referencia
a
la
librería
(http://docs.scipy.org/doc/scipy/reference/spatial.distance.html)
In [89]: from scipy.spatial.distance import pdist, squareform from scipy.cluster.hierarchy import linkage, dendrogram In [90]: np.random.seed(12345) clase1 = np.random.normal(loc=1,scale=0.2,size=(4,2)) clase2 = np.random.normal(loc=2,scale=0.2,size=(4,2)) clase3 = np.random.normal(loc=3,scale=0.2,size=(4,2)) clases = vstack((clase1, clase2, clase3)) x = clases[:,0] y = clases[:,1] plt.scatter(x,y, s=60) for i in range(12): plt.text(x[i]-0.025, y[i]+0.1,i)
Vamos a crear un dataframe con las 12 observaciones:
scipy.spatial.distance
In [109]: df = pd.DataFrame(clases, columns=['x', 'y']) df Out[109]:
x
y
0
0.959058 1.095789
1
0.896112 0.888854
2
1.393156 1.278681
3
1.018582 1.056349
4
2.153805 2.249287
5
2.201438 1.740756
6
2.054998 2.045783
7
2.270583 2.177286
8
2.599673 2.925631
9
3.333805 2.912286
10 2.892052 3.095397 11 3.649789 2.795754
La función pdist() calcula la distancia de cada uno de los 12 puntos con respecto a los demás. Se crea un array de
valores
In [93]: dm = pdist(df,metric='euclidean') dm.shape Out[93]: (66,)
A continuación, hacemos el clustering
In [101]: z = linkage(dm, method='complete') z.shape Out[101]: (11, 4) Y construimos un dendrograma
In [102]: dendrogram(z);
In [106]: # Se pueden buscar otras orientaciones dendrogram(z, orientation='right');
Por curiosidad, las distancias también se pueden poner en forma de matriz cuadrada, aunque como se ha visto, no es preciso para calcular el clustering:
In [94]: distxy = squareform(dm) distxy.shape Out[94]: (12, 12) In [95]: # Como vemos, es una matriz simétrica con 0 en la diagonal distxy[0:3,0:3] Out[95]: array([[ 0. , [ 0.21629657, [ 0.47105247,
k-clustering
0.21629657, 0. , 0.63167861,
0.47105247], 0.63167861], 0. ]])
El algoritmo "k-means" toma como entrada el número de clusters a generar, k (esto es su principal limitación) y un conjunto de vectores resultado de observaciones (en nuestro caso los 12 pares de coordenadas x,y). Devuelve un conjunto de k centroides, uno por cada cluster. Las observaciones son clasificadas mediante el número del cluster (el index del centroide más próximo). Este proceso se conoce a veces como "cuantificación" de los vectores de las observaciones. Al cluster index de un vector se le llama el "código" y la tabla que asocia códigos y centroides se conoce como el "code book"
In [107]: from scipy.cluster.vq import kmeans, vq In [127]: codebook, varianza = kmeans(df,3) -------------------------------------------------------------------------TypeError Traceback (most recent call la st) in () ----> 1 codebook, varianza = kmeans(df,3) /usr/lib/python2.7/dist-packages/scipy/cluster/vq.pyc in kmeans(obs, k_o r_guess, iter, thresh) 505 for i in range(iter): 506 #the intial code book is randomly selected from obse rvations --> 507 guess = take(obs, randint(0, No, k), 0) 508 book, dist = _kmeans(obs, guess, thresh = thresh) 509 if dist < best_dist: /usr/local/lib/python2.7/dist-packages/numpy/core/fromnumeric.pyc in tak e(a, indices, axis, out, mode) 116 except AttributeError: 117 return _wrapit(a, 'take', indices, axis, out, mode) --> 118 return take(indices, axis, out, mode) 119 120 TypeError: take() takes at most 4 arguments (5 given) In [128]: # parece que a la función kmeans no le gustan los dtaframes de pandas codebook, varianza = kmeans(np.vstack(zip(x,y)),3) In [129]: # sin embargo... pasandolo a array de numpy, funciona codebook, varianza = kmeans(df.values,3) In [132]: # En realidad, esto no lo necesitamos para nada print distortion 0.262435264634
In [131]: # estos son los centroides de los tres grupos: print codebook [[ 1.0667271 [ 2.17020602
1.07991825] 2.05327779]
[ 3.11882952
2.93226726]]
In [133]: # A continuación, la función vq() asigna números de clusters (códigos d el codebook) a las observaciones: In [134]: code,distance = vq(df.values,codebook) In [136]: # Así obtenemos el código de cada observación code Out[136]: array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]) In [138]: colores = ['red', 'blue', 'green'] c = [colores[i] for i in code] print c ['red', 'red', 'red', 'red', 'blue', 'blue', 'blue', 'blue', 'green', 'g reen', 'green', 'green']
In [139]: plt.scatter(x,y, s=60, c=c) Out[139]:
Ejemplo completo de clustering Para este ejemplo tenemos que leer un fichero de datos en el formato binario propietario de R: RData. Por eso lo mejor es leerlo con R y despues pasarlo a un dataframe de Pandas
In [148]: !ls ./datos foods-2011-10-03.json hola phoenix.html phoenix-tidied.html
prueba prueba2.csv prueba3.csv prueba4.csv
prueba4.txt prueba5.txt prueba.csv prueba.html
samsungData.rda warnings
In [144]: %load_ext rmagic In [149]: %%R load("./datos/samsungData.rda") write.csv(samsungData,file="./datos/samsungData.csv") In [150]: !ls ./datos foods-2011-10-03.json hola phoenix.html phoenix-tidied.html
prueba prueba2.csv prueba3.csv prueba4.csv
prueba4.txt prueba5.txt prueba.csv prueba.html
samsungData.csv samsungData.rda warnings
In [151]: samsungData = pd.read_csv('./datos/samsungData.csv') In [152]: samsungData.tail() Out[152]: Int64Index: 5 entries, 7347 to 7351 Columns: 564 entries, Unnamed: 0 to activity dtypes: float64(561), int64(2), object(1) In [153]: samsungData.shape Out[153]: (7352, 564) In [165]: # Nombres de las 10 primeras columnas samsungData.columns[0:10] Out[165]: Index([u'tBodyAcc-mean()-X', u'tBodyAcc-mean()-Y', u'tBodyAcc-mean()-Z', u'tBodyAcc-std()-X', u'tBodyAcc-std()-Y', u'tBodyAcc-std()-Z', u'tBodyA cc-mad()-X', u'tBodyAcc-mad()-Y', u'tBodyAcc-mad()-Z', u'tBodyAcc-max()X'], dtype=object) In [168]: # Nombres de las 10 últimas columnas samsungData.columns[-10:] Out[168]: Index([u'fBodyBodyGyroJerkMag-kurtosis()', u'angle(tBodyAccMean,gravity) ', u'angle(tBodyAccJerkMean),gravityMean)', u'angle(tBodyGyroMean,gravit yMean)', u'angle(tBodyGyroJerkMean,gravityMean)', u'angle(X,gravityMean) ', u'angle(Y,gravityMean)', u'angle(Z,gravityMean)', u'subject', u'activ ity'], dtype=object) In [161]: samsungData = samsungData.drop('Unnamed: 0', axis=1)
In [162]: samsungData.columns[0:10] Out[162]: Index([u'tBodyAcc-mean()-X', u'tBodyAcc-mean()-Y', u'tBodyAcc-mean()-Z', u'tBodyAcc-std()-X', u'tBodyAcc-std()-Y', u'tBodyAcc-std()-Z', u'tBodyA cc-mad()-X', u'tBodyAcc-mad()-Y', u'tBodyAcc-mad()-Z', u'tBodyAcc-max()X'], dtype=object) In [163]: samsungData['activity'].value_counts() Out[163]: laying 1407 standing 1374 sitting 1286 walk 1226 walkup 1073 walkdown 986 dtype: int64 In [178]: # para hacernos una idea de como viene codificado el dataframe # vamos a listar un subconjunto: samsungData.ix[985:995,[0,1,2,3,4,5,-2,-1]] Out[178]:
tBodyAcc- tBodyAcc- tBodyAcc- tBodyAcc- tBodyAcc- tBodyAcc subject mean()-X mean()-Y mean()-Z std()-X std()-Y std()-Z 985 0.198992
-0.002455
-0.117281
0.100775
0.355080
-0.266647
5
986 0.158101
-0.040474
-0.134750
0.063741
0.265524
-0.299606
5
987 0.281287
-0.034803
-0.089352
-0.064575
0.314012
-0.280909
5
988 0.418010
-0.016577
-0.153921
-0.047141
0.283657
-0.178543
5
989 0.428925
-0.037568
-0.169470
-0.033747
0.300829
-0.229894
5
990 0.292996
-0.036746
-0.111782
-0.953571
-0.863929
-0.870786
6
991 0.276552
-0.028512
-0.110449
-0.987560
-0.945003
-0.944290
6
992 0.271818
-0.032274
-0.113994
-0.995988
-0.959353
-0.955563
6
993 0.275229
-0.010966
-0.089999
-0.995814
-0.958768
-0.976571
6
994 0.279222
-0.005795
-0.092436
-0.996173
-0.969167
-0.980864
6
995 0.276892
-0.018711
-0.109727
-0.994897
-0.972814
-0.963744
6
Hay una línea por sujeto y lectura de los accelerómetros, y a cada vector de observaciones se asigna una actividad. A continuación vamos a ir probando variables del sujeto 1, y viendo si estas variables discriminan bien entre actividades:
In [205]: # Array de actividades acts = samsungData['activity'].unique() # secuencia de colores cols = 'bgrcmy'
In [206]: # Crear un diccionario de colores dic_col = {acts[i]:cols[i] for i in range(len(acts))} dic_col Out[206]: {'laying': 'r', 'sitting': 'g', 'standing': 'b', 'walk': 'c', 'walkdown': 'm', 'walkup': 'y'} In [207]: # Seleccionamos las filas del primer sujeto subj1 = samsungData[samsungData['subject']==1] # Ahora creamos un objeto "groupby" para agrupar por actividad grouped = subj1.groupby('activity') In [212]: #En abcisas vamos a representar los valores de la primera variable # Y en ordenadas el número de la observación (la lectura) # Siempre referido al sujeto 1 fig, (ax1, ax2) = plt.subplots(1,2, sharey=True) fig.set_size_inches(10, 5)
for act, df in grouped: ax1.scatter(df.ix[:,0], df.index, c=dic_col[act], label=act) ax2.scatter(df.ix[:,1], df.index, c=dic_col[act], label=act) ax1.set_xlabel(samsungData.columns[0]) ax1.set_ylabel(u'# Observación') ax2.set_xlabel(samsungData.columns[1]) ax2.legend(loc='upper left')
Conclusión: las dos primeras variables no nos permiten discriminar por tipos de actividad
Ahora vamos a probar a hacer un agrupamiento jerárquico basado en las tres primeras columnas, a ver si esto nos permite separar por grupos con un tipo de actividad
In [233]: # Construimos una variable categórica (factor) con las actividades # Ya que de esta manera actF.labels contendrá la actividad # codificada numéricamente de 0 a 5 actF = pd.Categorical.from_array(subj1['activity']) In [243]: dm = pdist(subj1.ix[:,0:3],metric='euclidean') z = linkage(dm, method='complete') # plt.figure(figsize=(5, 10)) # Una forma de dar el tamaño de la fi gura dendrogram(z, orientation='right', color_threshold=0.2, leaf_label_func=lambda n : 'X' * (actF.labels[n] + 1)); fig = plt.gcf() # Get current figure fig.set_size_inches(5,10) # Otra forma de dar el tamaño
Vamos a probar ahora con las columnas 9 y 10
