In [2]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
%matplotlib inline
On importe les données correspondant au réseau de vélos partagés lyonnais, et on considère cinq semaines de mesures, soit entre le lundi 21 août 2017 et le dimanche 25 septembre 2017.
In [3]:
lyon = pd.read_csv("velodata/lyon.csv", parse_dates = ['last_update'])
lyon = lyon.query("last_update >= '2017-08-21' and last_update < '2017-09-25'")
In [4]:
print(lyon.shape)
lyon.head(10)
Out[4]:
Pour se représenter l'activité d'une station Velov, on affiche une simple chronique temporelle du nombre de vélo disponible sur la station, ainsi que du nombre de stations disponibles. En principe, ces deux grandeurs donnent le nombre d'emplacements total sur la station, mais il peut arriver qu'un ou plusieurs emplacements soient en maintenance ; on conserve tout de même l'idée que le nombre total d'emplacement est une borne supérieure de la somme des vélos et des emplacements disponibles.
In [5]:
def plot_station(data, index):
"""Plot available bikes and bike stands for a given station"""
station = data.query("number == @index")
f, ax = plt.subplots(2, 1, figsize=(16, 6))
ax[0].plot(station.last_update, station.available_bike_stands, 'r-')
ax[0].set_xlabel("Available bike stands")
ax[1].plot(station.last_update, station.available_bikes, 'b-')
ax[1].set_xlabel("Available bikes")
In [6]:
plot_station(lyon, 3085)
Pour la prochaine section (analyse de séries temporelles), on ne considère qu'une station Velov, mais la même analyse peut être réalisée pour n'importe quelle station.
In [7]:
station = lyon.query("number == 3085")
avl_velos = station.set_index("last_update")[['available_bikes', 'bike_stands']].resample("15T")
avl_stands = station.set_index("last_update")[['available_bike_stands', 'bike_stands']].resample("15T")
Une première tentative d'analyse des séries temporelles est réalisée avec le package pyflux, solution standard dans cette optique.
In [8]:
import pyflux as pf
In [26]:
model = pf.ARIMA(avl_velos.reset_index(), 4, 4)
In [27]:
x = model.fit()
x.summary()
On obtient ce résultat, donnant la valeur des coefficients propres au modèle ARIMA. Trouver un jeu de paramètres pour lequel la matrice Hessienne est inversible reste toutefois compliqué, et dans ce cas de figure, il est impossible de confronter les valeurs obtenues à des tests statistiques.
On observe également que les coefficients associés à la partie auto-regressive du modèle sont nuls.
In [28]:
from fbprophet import Prophet
In [29]:
avl_velos = avl_velos.reset_index()
Inconvénient, il faut absolument que la colonne correspondant aux dates soit nommée ds, et que celle référant aux données mesurées soit nommée y. Une borne supérieure correspondant à la capacité du phénomène étudiée peut être rajoutée, la colonne correspondante devant être nommée cap. Ici, nous sommes dans ce cas, puisque la quantité de vélos disponible par station est bornée par le nombre d'emplacements.
In [30]:
avl_velos.columns = ["ds", "y", "cap"]
data = avl_velos.query("ds < '2017-09-18'")
future = avl_velos.query("ds >= '2017-09-18'")
Un modèle est déclaré selon un mode de croissance linear ou logistic, les différentes composantes de saisonnalité peuvent être précisées également (sans compter que des éléments de saisonnalité supplémentaires peuvent être rajoutés par la suite, cf API).
In [31]:
m = Prophet(growth='logistic',
daily_seasonality=True, weekly_seasonality=True, yearly_seasonality=False,
n_changepoints=10)
m.fit(data)
forecast = m.predict(future)
In [54]:
forecast.sample().T
Out[54]:
Prophet dispose d'une fonction d'affichage très propre, que l'on peut hacker via l'API de matplotlib (exemple ici : on rajoute les points correspondants aux observations réelles, en vert, à comparer aux prévisions du modèle, en bleu).
In [32]:
mp = m.plot(forecast)
test = mp.get_axes()
test[0].plot(future.ds, future.y, 'go', ms=3)
Out[32]:
Une autre fonction bien pratique permet d'afficher les composantes de la prévision, en séparant la courbe tendancielle des effets de saisonnalité.
In [53]:
m.plot_components(forecast)
Out[53]:
Malheureusement (au-delà du petit bug qui fait afficher les courbes en double), aucune fonction ne permet de récupérer les données ayant conduit à la génération de ces graphiques... Celle-ci peuvent malgré tout être reconstruites à la main, ce que nous pouvons montrer dans la section suivante.
Deux effets de saisonnalité sont cruciaux ici : l'effet de heure de la journée, et celui du jour de la semaine.
In [40]:
days_of_week = ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"]
forecast['day_index'] = forecast['ds'].dt.dayofweek
forecast['day'] = [days_of_week[index] for index in forecast['day_index']]
In [41]:
forecast.groupby("day")['weekly'].describe()
Out[41]:
In [44]:
f, ax = plt.subplots(1, 2, figsize=(16, 6))
sns.boxplot(x=forecast.day, y=forecast.weekly, ax=ax[0])
ax[1].plot(forecast.day_index.drop_duplicates(), forecast.groupby("day_index")["weekly"].mean(), '-', color='#0072B2')
f.tight_layout()
Le résultat sur les jours de semaine montre une utilisation plus intensive des Velov en semaine, et plus particulièrement le lundi. On note toutefois que ce résultat n'est pas équivalent à celui proposé via l'utilisation des fonctions d'affichage du package Prophet (alors même que le calcul à la main utilise la prédiction réalisée via Prophet)...
In [45]:
forecast['hour'] = forecast['ds'].dt.hour
In [46]:
forecast.groupby("hour")['daily'].describe()
Out[46]:
In [47]:
f, ax = plt.subplots(1, 2, figsize=(16, 6))
sns.boxplot(x=forecast.hour, y=forecast.daily, ax=ax[0])
ax[1].plot(forecast.hour.drop_duplicates(), forecast.groupby("hour")["daily"].mean(), '-', color='#0072B2')
f.tight_layout()
Le résultat en fonction des heures est aussi extrêmement intéressant : les Velov sont très utilisés en journée, et particulièrement sur des pointes matinales à 10h le matin et 15h l'après-midi. La nuit, la station tend généralement à être remplie. Au contraire des jours de la semaine, ce résultat correspond à celui proposé par le package Prophet.
On cherche à reproduire le travail précédent sur l'ensemble des stations pour caractériser chacune avec des profils temporels.
In [76]:
station_ids = lyon.number.drop_duplicates()
In [73]:
def daily_pattern(data, station):
avl_velos = data.query("number==@station").set_index("last_update")['available_bikes'].resample("15T").reset_index()
avl_velos['hour'] = avl_velos['last_update'].dt.hour
daily_pattern = avl_velos.groupby("hour")["available_bikes"].mean()
return daily_pattern
In [75]:
def weekly_pattern(data, station):
avl_velos = data.query("number==@station").set_index("last_update")['available_bikes'].resample("15T").reset_index()
avl_velos['day'] = avl_velos['last_update'].dt.dayofweek
weekly_pattern = avl_velos.groupby("day")["available_bikes"].mean()
return weekly_pattern
In [78]:
daily_patterns = pd.DataFrame([daily_pattern(lyon, station_id) for station_id in station_ids])
daily_patterns.index = station_ids
weekly_patterns = pd.DataFrame([weekly_pattern(lyon, station_id) for station_id in station_ids])
weekly_patterns.index = station_ids
In [86]:
from sklearn.cluster import KMeans
In [88]:
model = KMeans(n_clusters=4, random_state=0)
In [89]:
result = model.fit(weekly_patterns)
In [108]:
weekly_pattern_labels = pd.DataFrame({"label": result.labels_}, index=weekly_patterns.index)
In [93]:
import folium
In [95]:
position = [45.750000, 4.850000] # Lyon lat, lon coordinates
In [145]:
mp = folium.Map(location=position, zoom_start=13, tiles='cartodbpositron')
In [98]:
locations = pd.read_csv("./data/lyon-stations.csv")
In [109]:
locations = locations.merge(weekly_pattern_labels, right_index=True, left_on='idstation')
Out[109]:
In [130]:
locations.query("idstation==3085")
Out[130]:
In [100]:
colors = sns.color_palette('Set1', 4)
hex_colors = colors.as_hex()
In [143]:
locations["color"] = [hex_colors[label] for label in locations.label]
for _, row in locations.iterrows():
plt.plot(row.x, row.y, 'o', color=row.color)
In [146]:
for _,row in locations.iterrows():
folium.CircleMarker(
location=[row['x'], row['y']],
radius=5,
popup=row['nom'],
color=hex_colors[row['label']],
fill=True,
fill_opacity=0.7,
fill_color=hex_colors[row['label']]
).add_to(mp)
In [147]:
mp
Out[147]: