Estadistica Con Python I 1v444y

Tutorial sobre estadística en Python: I Distribuciones, gráficos, clustering In [1]: import pandas as pd import numpy as np import matplotlib.pyplot as plt import scipy.stats as s In [2]: print np.__version__ import scipy print scipy.__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:  (http://www.randalolson.com/2012/08/06/statistical­analysis­made­easy­in­python/) Esta  es  la  documentación  de  la  librería  (http://docs.scipy.org/doc/scipy/reference/tutorial/stats.html)

estadística 

de 

tutorial 

SciPy: 

R. 

Olson

scipy.stats

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  ree (https://github.com/herrfz/dataanalysis/blob/master/REE.md) El  material  completo  del  curso  de  Cursera  está  aquí:  Coursera  Data  Analysis (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]: # generación de valores aleatorios según distribución normal # normal(loc=0.0, scale=1.0, size=None) # loc es la media, scale la desviación estándar # size es un un valor entero # o una tupla con el shape deseado alturas = np.random.normal(188, 3, 20) print alturas [ 190.09944441 187.08032995 196.58267576 191.81629275

193.31053145 190.53477337 187.70675913 185.77887953

192.10140399 188.38027873 191.99139412 189.25968599

189.08948296 190.74397206 188.14624532 183.90044742

194.01804315 191.10964306 190.57930883 187.989887 ]

In [4]: # Cuando deseamos valores según la distribución normal estándar # (mean=0, stdev =1), podemos utilizar # standard_normal(size), donde size puede ser un valor entero o un shap e np.random.standard_normal((4,2)) Out[4]: array([[-0.12656278, 1.04534887], [ 0.96472613, -0.95483344], [ 0.40339144, -2.29152082], [-1.15241441, 0.73026361]]) In [22]: # Generación de valores aleatorios según distribución binomial # binomial(n, p, size) coinFlips = np.random.binomial(1, 0.5, 10) print coinFlips [0 1 1 1 1 1 0 0 1 1]

In [9]: # Generación de valores aleatorios según la distribución uniforme # En un intervalo semiabierto [low, high) # uniform(low=0.0, high=1.0, size=1) # size puede ser un entero o una tupla indicando un shape np.random.uniform(size=10) Out[9]: array([ 0.5526518 , 0.59462721, 0.13537618, 0.91364061, 0.59570865, 0.24854146, 0.7054177 , 0.17873533, 0.70354042, 0.30468929])

In [15]: # Generación de valores aleatorios enteros según la distribución unifor me # en el intervalo semiabierto [low, high) # o bien, si high se omite, en [0, low) # size es un entero o una tupla indicando un shape #randint(low, high=None, size=None) print np.random.randint(0,12,(3,4)) # También se puede utilizar random_integers, similar a la anterior # pero que genera valores enteros en el intervalo cerrado #[low, high] print print np.random.random_integers(1,6, size=10) [[ 5 2 0 8] [ 8 11 11 7] [11 9 0 9]] [6 6 3 4 1 4 1 6 6 1]

Obtención de muestras con Numpy In [16]: # random.choice(a, size=None, replace=True, p=None) # a es un array 1-D obligatoriamente, o bien un num. entero # (ver ejemplo siguiente) # size es un valor entero o una tupla indicando un shape # p es un array de probabilidades de igual longitud que a # si se omite se supone distribución uniforme para a # replace=False si se desea una muestra sin remplazamiento np.random.choice(alturas, size=5, replace=False) Out[16]: array([ 196.58267576, 190.74397206, 189.25968599, 191.99139412, 189.08948296]) In [18]: # Otra posibilidad es hacer a igual a un valor entero n # Viene bien para por ejemplo indexar con él # En este caso las muestras se toman del array arange(n) 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 5 1 4] [ 9 5 2 1] [ 6 1 11 9]] media total: 4.66666666667 media calculada variando la fila: [ 5.66666667 3.66666667 4.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 5.33333333 30.33333333 16.33333333] [ 3.33333333 12.91666667 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 5.50757055 4.04145188] [ 1.82574186 3.59397644 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

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  r,  y  se  interpreta  como  el porcentaje de variabilidad de la variable y  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  X ,  el  rango  de  X ,  RX   es  el  conjunto  de  valores  que  puede  tomar.  Cada valor  x ∈ RX   es  llamado  un  cuantil,  y  la  probabilidad  de  que  X   tome  el  valor  x   se  designa  por  P (X = x)

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, 5.71139329, 5.58134749, 6.19376692, 5.92728227, 4.37584035, 4.60381517, 7.66738297, 8.44622194, 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, 5.48459398, 3.25193664, 3.58073744, 7.89030591, 6.35598259, 4.12032957, 4.99845849, 6.55752111, 3.26517431])

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 , 0.25])

Gráficos exploratorios Representar muchos puntos

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]: <matplotlib.collections.PathCollection at 0x268a7950>

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. 0. 7. [ 0. 4. 26. [ 0. 10. 62. [ 2. 17. 134. [ 4. 19. 130. [ 0. 9. 93. [ 1. 4. 24. [ 1. 0. 6. [ 0. 0. 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. 0. 10. 6. 61. 16. 220. 67. 477. 135. 447. 147. 253. 72. 78. 20. 15. 6. 2. 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. -3.2 -2.4 -1.6 -0.8 0. [-4. -3.2 -2.4 -1.6 -0.8 0.

0.8 1.6 2.4 3.2 4. ] 0.8 1.6 2.4 3.2 4. ]

Ahora vamos a ver como podemos representar 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. 115.] [ 100. 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.1617511 ] [ 0.14065313 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

scipy.spatial.distance

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:

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 n(n − 1)/2  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], [ 0.21629657, 0. , 0.63167861], [ 0.47105247, 0.63167861, 0. ]])

k­clustering

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 <module>() ----> 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 1.07991825] [ 2.17020602 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]: <matplotlib.collections.PathCollection at 0x5927650>

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 samsungData.rda prueba5.txt warnings prueba.csv prueba.html

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 samsungData.csv prueba5.txt samsungData.rda prueba.csv warnings prueba.html

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 activ 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

walku

986 0.158101

­0.040474

­0.134750

0.063741

0.265524

­0.299606

5

walku

987 0.281287

­0.034803

­0.089352

­0.064575

0.314012

­0.280909

5

walku

988 0.418010

­0.016577

­0.153921

­0.047141

0.283657

­0.178543

5

walku

989 0.428925

­0.037568

­0.169470

­0.033747

0.300829

­0.229894

5

walku

990 0.292996

­0.036746

­0.111782

­0.953571

­0.863929

­0.870786

6

stand

991 0.276552

­0.028512

­0.110449

­0.987560

­0.945003

­0.944290

6

stand

992 0.271818

­0.032274

­0.113994

­0.995988

­0.959353

­0.955563

6

stand

993 0.275229

­0.010966

­0.089999

­0.995814

­0.958768

­0.976571

6

stand

994 0.279222

­0.005795

­0.092436

­0.996173

­0.969167

­0.980864

6

stand

995 0.276892

­0.018711

­0.109727

­0.994897

­0.972814

­0.963744

6

stand

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

In [245]: #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[:,9], df.index, c=dic_col[act], label=act) ax2.scatter(df.ix[:,10], df.index, c=dic_col[act], label=act) ax1.set_xlabel(samsungData.columns[9]) ax1.set_ylabel(u'# Observación') ax2.set_xlabel(samsungData.columns[10]) ax2.legend(loc='upper left')

Continuaremos este ejemplo más adelante una vez hayamos visto el análisis de componentes principales (PCA) In []: