Use a supervised machine learning algorithm to predict the availability for each bicyle-sharing stations in Lyon (France) based on the history data.
I use the tree method XGBoost to predict a "probability" of bikes availability for each station. A number close to 1. means that you have several available bikes. A number close to 0. means you don't have many bikes.
In [1]:
%matplotlib inline
In [32]:
import numpy as np
import pandas as pd
In [157]:
from xgboost import plot_tree
In [93]:
import matplotlib as mpl
from matplotlib import pyplot as plt
In [3]:
import folium
In [4]:
%load_ext watermark
In [25]:
%watermark -d -v -p numpy,pandas,xgboost,matplotlib,folium -g -m -w
The module prediction.py contains some functions dedicated to the bicyle-sharing stations predictions.
In [7]:
from prediction import (datareader, complete_data, cleanup, bikes_probability,
time_resampling, prepare_data_for_training, fit, prediction)
In [8]:
DATAFILE = './data/lyon.csv'
In [9]:
raw = datareader(DATAFILE)
In [10]:
raw.head()
Out[10]:
Min and max dates of the timeseries
In [11]:
print(raw.last_update.min())
print(raw.last_update.max())
Clean up some columns : drop some lines with the 'CLOSED' status, drop duplicates, remove missing values, etc.
In [12]:
df_clean = cleanup(raw)
Pipe some processing data functions :
10T)num_avail_bikes / total
In [13]:
df = (df_clean.pipe(time_resampling)
.pipe(complete_data)
.pipe(bikes_probability))
In [14]:
df.head()
Out[14]:
This is the final dataset. For further prediction, I could add some weather forecasts data to these features.
Let's select a time window (start, stop) to a single prediction.
In [21]:
start = pd.Timestamp("2017-07-11T00:00:00") # Tuesday
predict_date = pd.Timestamp("2017-07-26T10:00:00") # wednesday
# predict the next 30 minutes
freq = '30T'
# number of predictions at 'predict_date'.
# Here, the next 30 minutes and the next hour (30 minutes + 30 minutes).
# If you want to predict the next 3 hours, every 30 minutes, thus set periods=6
periods = 2
In [22]:
train_X, train_Y, test_X, test_Y = prepare_data_for_training(df,
predict_date,
freq=freq,
start=start,
periods=periods,
observation='probability')
In [23]:
train_X.head()
Out[23]:
In [24]:
# Observation with a shift of T+30 minutes in comparison to train_X.
# This is the 'future' availability used for traning.
train_Y.head()
Out[24]:
In [28]:
train_X.iloc[3200:3210]
Out[28]:
The fit function create some data structure for the XGBoost from the train and test DataFrames (i.e. xgb.DMatrix), configure the model and launch it with the objective: 'reg:logistic'. It's a regression, not a classification.
In [30]:
# the 'booster'
bst = fit(train_X, train_Y, test_X, test_Y)
In [31]:
# compute the prediction from test_*
pred = prediction(bst, test_X, test_Y)
In [34]:
pred[:5]
Out[34]:
In [37]:
print("Number of predictions: {}".format(len(pred)))
In [33]:
# Compute the RMSE
rmse = np.sqrt(np.mean((pred - test_Y)**2))
rmse
Out[33]:
In [159]:
# must install graphviz
# plot_tree(bst)
In [35]:
result = test_X.copy()
result['ts_future'] = test_Y.index.shift(1, freq=freq)
result['observation'] = test_Y.copy()
result['ts_future'] = test_Y.index.shift(1, freq=freq)
result['prediction'] = pred
result['error'] = pred - test_Y
result['relative_error'] = 100. * np.abs(pred - test_Y) / test_Y
result['quad_error'] = (pred - test_Y)**2
result.to_csv("prediction-freq-{}-{}.csv".format(freq, predict_date))
In [36]:
result.head(10)
Out[36]:
CSV file with station coordinates
In [68]:
locations = pd.read_csv("./data/lyon-stations.csv")
In [69]:
locations.shape
Out[69]:
Some stations were removed when the data were cleaned up. Remove them from the location data.
In [73]:
mask = locations['idstation'].isin(result.station.unique())
In [74]:
mask.sum()
Out[74]:
In [75]:
locations = locations[mask]
In [76]:
locations = locations.rename_axis({'idstation': 'station'}, axis=1)
In [77]:
locations.head()
Out[77]:
Some station names contains the ' character. Replace it by the HTML code for folium.
In [78]:
locations["nom"] = locations['nom'].str.replace("'", "'")
Select the prediction data for a specific timestamp
In [79]:
data_to_plot = result.loc[predict_date]
In [80]:
data_to_plot.shape
Out[80]:
In [81]:
data_to_plot.head()
Out[81]:
In [83]:
yhat = data_to_plot[['station', 'prediction']].merge(locations, on='station')
yhat.head()
Out[83]:
In [104]:
y = data_to_plot[['station', 'observation']].merge(locations, on='station')
In [140]:
error = data_to_plot[['station', 'error']].merge(locations, on='station')
In [122]:
colormap = 'RdYlBu'
cmap = plt.get_cmap(colormap)
In [137]:
# show the colormap use to plot the stations, values [0, 1]
gradient = np.linspace(0, 1, 256)
gradient = np.vstack((gradient, gradient))
fig, ax = plt.subplots(1)
fig.subplots_adjust(top=0.95, bottom=0.80, left=0.2, right=0.99)
ax.set_xticks([0., 64, 128, 192, 256])
ax.set_xticklabels([0., 0.25, 0.5, 0.75, 1.])
ax.set_xlabel('Bikes Availability')
ax.imshow(gradient, aspect='auto', cmap=cmap, vmin=0, vmax=1)
plt.title('Colormap used to plot stations')
Out[137]:
In [95]:
color = lambda x: mpl.colors.to_hex(cmap(x))
In [89]:
# Lyon (France) Position
position = [45.750000, 4.850000]
In [106]:
mp_pred = folium.Map(location=position, zoom_start=13, tiles='cartodbpositron')
In [107]:
# Map of the predicted values
for _,row in yhat.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=2,
popup=row['nom'],
color=color(row['prediction']),
fill=True,
fill_opacity=0.3,
fill_color=color(row['prediction'])
).add_to(mp_pred)
In [102]:
mp_pred
Out[102]:
In [108]:
# Map for the observation
mp_obs = folium.Map(location=position, zoom_start=13, tiles='cartodbpositron')
In [109]:
# Map of the observations
for _,row in y.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=2,
popup=row['nom'],
color=color(row['observation']),
fill=True,
fill_opacity=0.3,
fill_color=color(row['observation'])
).add_to(mp_obs)
In [138]:
mp_obs
Out[138]:
In [153]:
# Colormap for error (by default, the color map fits for [0, 1] values)
norm = mpl.colors.Normalize(vmin=-1, vmax=1)
color_error = lambda x: mpl.colors.to_hex(cmap(norm(x)))
In [154]:
# Map for the errors
mp_error = folium.Map(location=position, zoom_start=13, tiles='cartodbpositron')
In [155]:
# Map of the errors
for _,row in error.iterrows():
folium.CircleMarker(
location=[row['lat'], row['lon']],
radius=2,
popup=row['nom'],
color=color_error(row['error']),
fill=True,
fill_opacity=0.3,
fill_color=color_error(row['error'])
).add_to(mp_error)
In [156]:
mp_error
Out[156]: