2A.ml - Clustering - correction

Links: notebook, html, PDF, python, slides, GitHub

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.

from jyquickhelper import add_notebook_menu
add_notebook_menu()
%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.

from pyensae.datasource import download_data
file = download_data("features_bike_chicago.zip")
file
['features_bike_chicago.txt']
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()
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

features.shape
(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.

names = features.columns[3:]
from sklearn.cluster import KMeans
clus = KMeans(10)
clus.fit(features[names])
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)
pred = clus.predict(features[names])
set(pred)
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
features["cluster"] = pred
features[["cluster", "weekday", "station_id"]].groupby(["cluster", "weekday"]).count()
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
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");
<matplotlib.axes._subplots.AxesSubplot at 0x149c1117278>
../_images/td2a_clustering_correction_14_1.png

Let’s draw the clusters.

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)
../_images/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

from sklearn.cluster import DBSCAN
dbs = DBSCAN(eps=0.1)
pred_dbs = dbs.fit_predict(features[names])
set(pred_dbs)
{-1, 0, 1, 2, 3, 4, 5}
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");
../_images/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.

piv = features.pivot_table(index=["station_id", "station_name"],
                           columns="weekday", values="cluster")
piv.head()
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
piv["distincts"] = piv.apply(lambda row: len(set(row[i] for i in range(0,7))), axis=1)
pivn = piv.reset_index()
pivn.head()
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
pivn.columns = [str(_).replace(".0", "") for _ in pivn.columns.values]
pivn.head()
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.

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")
stations = pandas.read_csv("Divvy_Stations_2016_Q3.csv")
stations.head()
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
data = stations.merge(pivn, left_on=["id", "name"],
                      right_on=["station_id", "station_name"], suffixes=('_s', '_c'))
data.sort_values("id").head()
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
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)
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%")
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