.. _td2aclusteringcorrectionrst: =============================== 2A.ml - Clustering - correction =============================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/td2a_ml/td2a_clustering_correction.ipynb|*` Ce notebook utilise les données des vélos de Chicago `Divvy Data `__. Il s’inspire du challenge créée pour découvrir les habitudes des habitantes de la ville `City Bike `__. L’idée est d’explorer plusieurs algorithmes de clustering et comment trafiquer les données pour les faire marcher. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: .. code:: ipython3 %matplotlib inline Les données ----------- Elles ont été prétraitées selon le notebook `Bike Pattern 2 `__. Elles représentent la distribution du nombre de vélos partant (*startdist*) et arrivant (*stopdist*). On utilise le clustering pour découvrir les différents usages des habitants de Chicago avec pour intuition le fait que les habitants de Chicago utilise majoritairement les vélos pour aller et venir entre leur appartement et leur lieu de travail. Cette même idée mais à Paris est illustrée par ce billet de blog : `Busy areas in Paris `__. .. code:: ipython3 from pyensae.datasource import download_data file = download_data("features_bike_chicago.zip") file .. parsed-literal:: ['features_bike_chicago.txt'] .. code:: ipython3 import pandas features = pandas.read_csv("features_bike_chicago.txt", sep="\t", encoding="utf-8", low_memory=False, header=[0,1]) features.columns = ["station_id", "station_name", "weekday"] + list(features.columns[3:]) features.head() .. raw:: html
station_id station_name weekday (startdist, 00:00:00) (startdist, 00:10:00) (startdist, 00:20:00) (startdist, 00:30:00) (startdist, 00:40:00) (startdist, 00:50:00) (startdist, 01:00:00) ... (stopdist, 22:20:00) (stopdist, 22:30:00) (stopdist, 22:40:00) (stopdist, 22:50:00) (stopdist, 23:00:00) (stopdist, 23:10:00) (stopdist, 23:20:00) (stopdist, 23:30:00) (stopdist, 23:40:00) (stopdist, 23:50:00)
0 2.0 Michigan Ave & Balbo Ave 0.0 0.003756 0.000939 0.002817 0.000000 0.000000 0.003756 0.000000 ... 0.004381 0.002191 0.004381 0.002191 0.004381 0.004381 0.005476 0.002191 0.000000 0.005476
1 2.0 Michigan Ave & Balbo Ave 1.0 0.000000 0.000000 0.001106 0.001106 0.001106 0.002212 0.000000 ... 0.009371 0.012048 0.006693 0.004016 0.005355 0.006693 0.002677 0.000000 0.000000 0.000000
2 2.0 Michigan Ave & Balbo Ave 2.0 0.001357 0.002714 0.000000 0.001357 0.000000 0.005427 0.000000 ... 0.002907 0.002907 0.015988 0.005814 0.001453 0.001453 0.011628 0.000000 0.000000 0.007267
3 2.0 Michigan Ave & Balbo Ave 3.0 0.000000 0.004144 0.000000 0.000000 0.002762 0.004144 0.000000 ... 0.009274 0.003091 0.003091 0.007728 0.001546 0.003091 0.009274 0.001546 0.007728 0.001546
4 2.0 Michigan Ave & Balbo Ave 4.0 0.000000 0.000000 0.000000 0.002846 0.000000 0.000000 0.000949 ... 0.008214 0.001027 0.006160 0.004107 0.015400 0.006160 0.002053 0.006160 0.007187 0.000000

5 rows × 291 columns

.. code:: ipython3 features.shape .. parsed-literal:: (4051, 291) Les données sont agrégrées par tranche de 10 minutes soit 144 période durant la journée et 288 nombre pour les départs et arrivées de vélos. Cela explique les dimensions de la matrice. k-means ------- On cherche à trouver différentes zones de la villes pour différents usages. Zones de résidences, zones de travail, zone d’amusements et on suppose que les heures de départs et d’arrivées reflètent ces usages. Dans une zone de travail, typiquement le quartier d’affaires, les vélos arriveront principalement le matin et repartiront le soir. Ce sera l’inverse pour une zone de résidences. C’est pour cela que les arrivées et les départs des vélos ont été agrégés par jour de la semaine. La distribution des arrivées risquent d’être bien différentes le week-end. .. code:: ipython3 names = features.columns[3:] .. code:: ipython3 from sklearn.cluster import KMeans clus = KMeans(10) clus.fit(features[names]) .. parsed-literal:: KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300, n_clusters=10, n_init=10, n_jobs=1, precompute_distances='auto', random_state=None, tol=0.0001, verbose=0) .. code:: ipython3 pred = clus.predict(features[names]) set(pred) .. parsed-literal:: {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} .. code:: ipython3 features["cluster"] = pred .. code:: ipython3 features[["cluster", "weekday", "station_id"]].groupby(["cluster", "weekday"]).count() .. raw:: html
station_id
cluster weekday
0 0.0 16
1.0 11
2.0 16
3.0 17
4.0 7
5.0 7
6.0 12
1 4.0 1
5.0 1
6.0 2
2 0.0 281
1.0 322
2.0 316
3.0 305
4.0 285
5.0 16
6.0 15
3 0.0 1
6.0 1
4 0.0 140
1.0 103
2.0 98
3.0 113
4.0 143
5.0 545
6.0 535
5 5.0 1
6 3.0 1
4.0 1
5.0 3
6.0 1
7 0.0 135
1.0 140
2.0 148
3.0 141
4.0 143
5.0 7
6.0 11
8 0.0 2
1.0 2
2.0 1
6.0 2
9 0.0 1
2.0 1
6.0 1
.. code:: ipython3 nb = features[["cluster", "weekday", "station_id"]].groupby(["cluster", "weekday"]).count() nb = nb.reset_index() nb[nb.cluster.isin([0, 3, 5, 6])].pivot("weekday","cluster", "station_id").plot(kind="bar"); .. parsed-literal:: .. image:: td2a_clustering_correction_14_1.png Let’s draw the clusters. .. code:: ipython3 centers = clus.cluster_centers_.T import matplotlib.pyplot as plt fig, ax = plt.subplots(centers.shape[1], 2, figsize=(10,10)) nbf = centers.shape[0] // 2 x = list(range(0,nbf)) col = 0 dec = 0 colors = ["red", "yellow", "gray", "green", "brown", "orange", "blue"] for i in range(centers.shape[1]): if 2*i == centers.shape[1]: col += 1 dec += centers.shape[1] color = colors[i%len(colors)] ax[2*i-dec, col].bar (x, centers[:nbf,i], width=1.0, color=color) ax[2*i-dec, col].set_ylabel("cluster %d - start" % i, color=color) ax[2*i+1-dec, col].bar (x, centers[nbf:,i], width=1.0, color=color) ax[2*i+1-dec, col].set_ylabel("cluster %d - stop" % i, color=color) .. image:: td2a_clustering_correction_16_0.png On a réussi à isoler plusieurs usages différents. On voit les départs les matin et les arrivées le soir, le modèle inversé, un autre… Mais les-a-t-on tous trouvés ? Exercice 1 : petits clusters ---------------------------- Les petits clusters viennent du fait qu’on impose un nombre de clusters à l’algorithme qui s’efforce de les trouver. Quelques points aberrents sont trop éloignés et isolés du reste des autres. Pour les enlever, on peut simplement réduire le nombre de clusters ou éliminer les points aberrants et les garder comme anomalies. Les petits clusters sont des stations peu utilisées. Le problème quand on calcule une distribution et qu’on la normalise, typiquement celles des arrivées de vélos à une stations, c’est qu’on perd l’information du nombre de vélos ayant servi à l’estimer. L’algorithme `k-means++ `__ commence à choisir les points les plus éloignés des uns des autres et cela tombe sur ces stations peu utilisées. Vu le nombre de dimension, 288, les distances sont rapidement très grandes. Une idée consiste à retirer ces stations peu utilisées ou boucler sur un algorithme du type : 1. Faire tourner les k-means. 2. Retirer les stations dans des clusters isolées. 3. Retour à l’étape 1 jusqu’à ce que cela soit interprétable. Exercice 2 : autres types de clustering --------------------------------------- .. code:: ipython3 from sklearn.cluster import DBSCAN dbs = DBSCAN(eps=0.1) pred_dbs = dbs.fit_predict(features[names]) set(pred_dbs) .. parsed-literal:: {-1, 0, 1, 2, 3, 4, 5} .. code:: ipython3 features["cluster_dbs"] = pred_dbs nbs = features[["cluster_dbs", "weekday", "station_id"]].groupby(["cluster_dbs", "weekday"]).count() nbs = nbs.reset_index() nbs[nbs.cluster_dbs.isin([0, 3, 5, 6])].pivot("weekday","cluster_dbs", "station_id").plot(kind="bar"); .. image:: td2a_clustering_correction_22_0.png L’algorithme `dbscan `__ utilise la proximité des points pour classer les observations. On se retrouve dans le cas où il n’y a pas vraiment de frontière entre les clusters et tous les points se retrouvent associés en un unique cluster excepté quelques points aberrants. L’algorithme `DBScan `__ ne fonctionne pas sur ces données. Une des raisons que la distance dans un espace vectoriels avec tant de dimensions n’est pas loin d’une information type binaire : différent partout ou identique presque partout. Qu’un vélo arrive 9h10 ou à 9h20, cela n’a pas beaucoup d’importance pour notre problème, pourtant, la distance choisie lui donnera autant d’importance que si le vélo était arrivé à 10h du soir. Pour éviter cela, il faudrait lisser les distributions avec une moyenne mobile. Graphes ------- Un peu de code pour voir la réparition des clusters sur une carte. .. code:: ipython3 piv = features.pivot_table(index=["station_id", "station_name"], columns="weekday", values="cluster") piv.head() .. raw:: html
weekday 0.0 1.0 2.0 3.0 4.0 5.0 6.0
station_id station_name
2.0 Michigan Ave & Balbo Ave 4.0 2.0 4.0 4.0 4.0 4.0 4.0
3.0 Shedd Aquarium 4.0 4.0 4.0 4.0 4.0 4.0 4.0
4.0 Burnham Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0
5.0 State St & Harrison St 2.0 2.0 2.0 2.0 4.0 4.0 4.0
6.0 Dusable Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0
.. code:: ipython3 piv["distincts"] = piv.apply(lambda row: len(set(row[i] for i in range(0,7))), axis=1) .. code:: ipython3 pivn = piv.reset_index() pivn.head() .. raw:: html
weekday station_id station_name 0.0 1.0 2.0 3.0 4.0 5.0 6.0 distincts
0 2.0 Michigan Ave & Balbo Ave 4.0 2.0 4.0 4.0 4.0 4.0 4.0 2
1 3.0 Shedd Aquarium 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
2 4.0 Burnham Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
3 5.0 State St & Harrison St 2.0 2.0 2.0 2.0 4.0 4.0 4.0 2
4 6.0 Dusable Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
.. code:: ipython3 pivn.columns = [str(_).replace(".0", "") for _ in pivn.columns.values] pivn.head() .. raw:: html
station_id station_name 0 1 2 3 4 5 6 distincts
0 2.0 Michigan Ave & Balbo Ave 4.0 2.0 4.0 4.0 4.0 4.0 4.0 2
1 3.0 Shedd Aquarium 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
2 4.0 Burnham Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
3 5.0 State St & Harrison St 2.0 2.0 2.0 2.0 4.0 4.0 4.0 2
4 6.0 Dusable Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
Une carte des stations un jour de semaine. .. code:: ipython3 from pyensae.datasource import download_data if False: # Provient du site de Chicago file = download_data("Divvy_Trips_2016_Q3Q4.zip", url="https://s3.amazonaws.com/divvy-data/tripdata/") else: # Copie au cas où celui-ci tomberait en panne file = download_data("Divvy_Trips_2016_Q3.zip") .. code:: ipython3 stations = pandas.read_csv("Divvy_Stations_2016_Q3.csv") stations.head() .. raw:: html
id name latitude longitude dpcapacity online_date
0 456 2112 W Peterson Ave 41.991178 -87.683593 15 5/12/2015
1 101 63rd St Beach 41.781016 -87.576120 23 4/20/2015
2 109 900 W Harrison St 41.874675 -87.650019 19 8/6/2013
3 21 Aberdeen St & Jackson Blvd 41.877726 -87.654787 15 6/21/2013
4 80 Aberdeen St & Monroe St 41.880420 -87.655599 19 6/26/2013
.. code:: ipython3 data = stations.merge(pivn, left_on=["id", "name"], right_on=["station_id", "station_name"], suffixes=('_s', '_c')) data.sort_values("id").head() .. raw:: html
id name latitude longitude dpcapacity online_date station_id station_name 0 1 2 3 4 5 6 distincts
357 2 Michigan Ave & Balbo Ave 41.872638 -87.623979 35 5/8/2015 2.0 Michigan Ave & Balbo Ave 4.0 2.0 4.0 4.0 4.0 4.0 4.0 2
456 3 Shedd Aquarium 41.867226 -87.615355 31 4/24/2015 3.0 Shedd Aquarium 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
53 4 Burnham Harbor 41.856268 -87.613348 23 5/16/2015 4.0 Burnham Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
497 5 State St & Harrison St 41.874053 -87.627716 23 6/18/2013 5.0 State St & Harrison St 2.0 2.0 2.0 2.0 4.0 4.0 4.0 2
188 6 Dusable Harbor 41.885042 -87.612795 31 4/24/2015 6.0 Dusable Harbor 4.0 4.0 4.0 4.0 4.0 4.0 4.0 1
.. code:: ipython3 def folium_html_stations_map(stations, html_width=None, html_height=None, **kwargs): import folium from pyensae.notebookhelper import folium_html_map map_osm = None for key, value in stations: x, y = key if map_osm is None: if "zoom_start" not in kwargs: kwargs["zoom_start"] = 11 if "location" not in kwargs: map_osm = folium.Map(location=[x, y], **kwargs) else: map_osm = folium.Map(kwargs["location"], **kwargs) if isinstance(value, tuple): name, value = value map_osm.add_child(folium.CircleMarker( [x, y], popup=name, radius=15, fill_color=value, color=value)) else: map_osm.add_child(folium.CircleMarker( [x, y], radius=15, fill_color=value, color=value)) return folium_html_map(map_osm, width=html_width, height=html_height) .. code:: ipython3 colors = ["red", "yellow", "gray", "green", "brown", "orange", "blue", "black", "pink", "violet"] for i, c in enumerate(colors): print("Cluster {0} is {1}".format(i, c)) xy = [] for els in data.apply(lambda row: (row["latitude"], row["longitude"], row["1"], row["name"]), axis=1): try: cl = int(els[2]) except: # NaN continue name = "%s c%d" % (els[3], cl) color = colors[cl] xy.append( ( (els[0], els[1]), (name, color))) folium_html_stations_map(xy, width="80%") .. parsed-literal:: Cluster 0 is red Cluster 1 is yellow Cluster 2 is gray Cluster 3 is green Cluster 4 is brown Cluster 5 is orange Cluster 6 is blue Cluster 7 is black Cluster 8 is pink Cluster 9 is violet .. raw:: html