План исследования
Подробности на kaggle.
from datetime import datetime as datetime
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, TimeSeriesSplit, cross_val_score, GridSearchCV
from sklearn.model_selection import validation_curve, learning_curve
from sklearn.metrics import r2_score, mean_absolute_error
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from xgboost import XGBRegressor
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from plotly import graph_objs as go
init_notebook_mode(connected = True)
import warnings
warnings.filterwarnings('ignore')
Интересующие нас данные хранятся в следующих файлах:
Целевой признак - количество поездок на великах в день
Загружаем данные из файлов
trip = pd.read_csv('trip.csv', parse_dates=['start_date', 'end_date'])
weather = pd.read_csv('weather.csv', parse_dates=['date'])
station = pd.read_csv('station.csv', parse_dates=['installation_date'])
Переведём продолжительность в минуты и посмотрим на распределение
trip['duration'] //= 60
trip[['duration']].describe(percentiles=[.10, .25, .50, .75, .90])
Ого... Максимальное значение почти 200 дней, хватит на пару кругостветок (если вы, например, этот парень). Думаю, что использование дольше нескольких дней можно смело отнести к выбросам, т.к. нас будут интересовать постоянные пользователи велопарковки.
Как правило люди берут велик, что бы доехать от одной парковки до другой, т.к. 90% пользуются им меньше 20 минут, чего (по моему сугубо личному мнению) маловато, чтобы насладиться поездкой, т.е. велосипеды используются как общественный транспорт, соответвенно спрос должен коррелировать с началом и окончанием рабочего дня и падать по выходным, позже мы проверим эту теорию.
Посмотрим как выглядят наши данные
trip.shape
trip.head()
Проверим пропуски
trip.isnull().sum()
У некоторых записей не заполнен почтовый код, но это не страшно, он нам не пригодится.
Дальше разберёмся с погодными условиями, посмотрим на распределения погодных условий, данных о ней много, даже слишком.
weather[[col for col in weather.columns if 'min' in col]].describe()
weather[[col for col in weather.columns if 'mean' in col]].describe()
weather[[col for col in weather.columns if 'max' in col]].describe()
weather.shape
Проверим пропуски
weather.isnull().sum()
С данными о погоде всё не так радужно, ноо... лучше, чем могло бы быть. Наибольшее количество пропусков в признаке events. Это можно объяснить отсутствием особенных погодных условий, то есть был обычный ясный день.
Сравним даты данных о погодных условиях и поездках
weather.date.min(), weather.date.max(), weather.date.max() - weather.date.min()
trip.start_date.min(), trip.end_date.max(), trip.end_date.max() - trip.start_date.min()
Промежутки совпадают, ничего лишнего нет, отлично.
weather.head()
Напоследок данные о пунктах выдачи
station.shape
trip.start_station_id.unique().shape[0] == station.shape[0]
Всего имеем 70 парковок и из каждой из них имеется информация в данных.
station.head()
Рассмотрим количество поездок за каждый день
trip_days = trip['start_date'].apply(lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()
trip_days.columns = ['trips']
def plotMovingAverage(df, n, title):
rolling_mean = df.rolling(window=n).mean()
trace0 = go.Scatter(
x = df.index,
y = df['trips'],
mode = 'lines',
name='Actual values'
)
trace1 = go.Scatter(
x = rolling_mean.index,
y = rolling_mean['trips'],
mode = 'lines',
name='Rolling mean trend'
)
data = [trace0, trace1]
layout = {'title': title}
fig = dict(data=data, layout=layout)
iplot(fig, show_link=False)
plotMovingAverage(trip_days, 7, 'Сглаживаем по неделям')
plotMovingAverage(trip_days, 30, 'Сглаживаем по месяцам')
Можно заметить сразу несколько интересных вещей:
Во-первых, что очевидно, поездки зависят от времени года, минимальное количество зимой, максимальное - летом.
Во-вторых, похоже наибольшей популярность велики пользуется в будние дни, и падает на выходных. Убедимся в этом построим гистограмму распределения по дням недели.
weekdays = {
0: 'Понедельник',
1: 'Вторник',
2: 'Среда',
3: 'Четверг',
4: 'Пятница',
5: 'Суббота',
6: 'Воскресенье'
}
trip_weekdays = trip['start_date'].apply(lambda dt: dt.weekday()).value_counts().sort_index().to_frame()
trip_weekdays.columns = ['trips']
trip_weekdays.index = trip_weekdays.index.map(lambda day: weekdays[day])
trace = go.Bar(
x = trip_weekdays.index,
y = trip_weekdays['trips']
)
data = [trace]
layout = {'title': 'Распределение по дням недели'}
fig = dict(data=data, layout=layout)
iplot(fig, show_link=False)
trip['is_weekend'] = trip['start_date'].apply(lambda dt: int(dt.weekday() >= 5))
trip_weekday = trip[trip['is_weekend'] == 0]['start_date'].apply(
lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()
trip_weekday.columns = ['trips']
trip_weekend = trip[trip['is_weekend'] == 1]['start_date'].apply(
lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()
trip_weekend.columns = ['trips']
trace0 = go.Scatter(
x = trip_weekday.index,
y = trip_weekday['trips'],
mode = 'markers',
name='Будни'
)
trace1 = go.Scatter(
x = trip_weekend.index,
y = trip_weekend['trips'],
mode = 'markers',
name='Выходные'
)
data = [trace0, trace1]
layout = {'title': 'Распределение по будням/выходным'}
fig = dict(data=data, layout=layout)
iplot(fig, show_link=False)
Теперь рассмотрим распределение поездок по часам
trip_hours = trip['start_date'].apply(lambda dt: dt.hour).value_counts().sort_index().to_frame()
trip_hours.columns = ['trips']
trace = go.Bar(
x = trip_hours.index,
y = trip_hours['trips']
)
data = [trace]
layout = {'title': 'Распределение по часам'}
fig = dict(data=data, layout=layout)
iplot(fig, show_link=False)
Итак по результатам предварительного анализа можно сделать следующие выводы:
Для начала посмотрим на количество выбросов и отфильтруем поездки более суток.
"Выбросы: {0:d} шт, {1:.2f}%".format((trip['duration'] >= 60 * 24).sum(),
(trip['duration'] >= 60 * 24).sum() / trip.shape[0] * 100)
trip = trip[trip['duration'] < 60 * 24]
Теперь начнём формирование обучающей выборки
# Просуммируем количество поездок по дням и отсортируем по дате
data = trip['start_date'].apply(lambda dt: datetime(dt.year, dt.month, dt.day, 0, 0)).value_counts().sort_index().to_frame()
# Приведём наименование
data.columns = ['target']
data['date'] = data.index
data.reset_index(drop=True, inplace=True)
#data = data.reindex(['date', 'target'], axis=1)
data.head()
Теперь добавим информацию о погодных условиях
Информация о погодных условиях разбита по почтовым кодам, так что очевидно данных о погоде больше, чем нужно
data.shape[0], weather.shape[0]
weather['zip_code'].unique()
Посмотрим на пропуски в погодных данных для разных кодов.
for code in weather['zip_code'].unique():
print("Код: {0}. Количество пропусков: {1}".format(code, weather[weather.zip_code == code].isnull().sum().values.sum()))
Данные с кодом 94107 имееют меньше всего пропусков. Велопарковки расположены относительно недалеко друг от друга, поэтому можно выбрать эти данные в качестве основных погодных условий без большого риска в потере качества предсказаний.
weather = weather[weather['zip_code'] == 94107]
Теперь разберёмся с пропусками в погодных данных
weather.isnull().sum()
Как уже говорилось выше, пропуск значения признака events означает отсутсвие особенных погодных условий. Заполним все пропуски этого признака новым погодным условием - Fair
.
weather.loc[weather.events.isnull(), 'events'] = "Fair"
weather['events'].unique()
Среди значений погодных условий есть дублирующие Rain
и rain
. Оставим только первый вариант простой заменой одного на другое.
weather.loc[weather['events'] == 'rain', 'events'] = "Rain"
weather.isnull().sum()
Отлично! Сам себя не похвалишь... Теперь уберём пропуски в признаке Максимальная скорость порыва ветра
.
Предположим, что скорость максимального порыва должна сильно коррелировать с максимальной скоростью ветра. Найдём медианные значения порыва ветра по скорости и заменим пропуски этим значением.
gust_by_wind = weather.groupby('max_wind_Speed_mph')['max_gust_speed_mph'].median()
def fill_gust_by_wind(row):
row['max_gust_speed_mph'] = gust_by_wind[row['max_wind_Speed_mph']]
return row
weather[weather['max_gust_speed_mph'].isnull()] = \
weather[weather['max_gust_speed_mph'].isnull()].apply(fill_gust_by_wind, axis=1)
weather.isnull().sum()
Посмотрим на типы данных в колонках информации о погоде
weather.dtypes
Признак precipitation_inches
имеет тип object, что немного противоречит природе вещей и крайне не нравится методам обучения, которые мы будем использовать. Исправим это и снова заполним получившиеся пропуски медианными значениями.
weather['precipitation_inches'] = pd.to_numeric(weather['precipitation_inches'], errors = 'coerce')
weather.loc[weather['precipitation_inches'].isnull(), 'precipitation_inches'] = \
weather[weather['precipitation_inches'].notnull()]['precipitation_inches'].median()
Остаётся только добавить информацию о велопарковках и собрать всё вместе в обучающую выборку.
station.head()
Посчитаем количество парковочных мест на каждый день из обучающей выборки
total_docks = []
for day in data['date']:
total_docks.append(sum(station[station.installation_date <= day].dock_count))
data['total_docks'] = total_docks
Теперь, когда мы закончили предварительную обработку данных сформируем признаки для обучающей выборки
weather = pd.concat([weather, pd.get_dummies(weather['events'])], axis=1)
# Выкидываем лишние признаки и добавляем в обучающую выборку
weather.drop(['events', 'zip_code'], axis=1, inplace=True)
data = data.merge(weather, left_on='date', right_on='date')
Так же добавим дополнительные временны́е признаки
data['weekend'] = data['date'].apply(lambda dt: int(dt.weekday() >= 5))
data['weekday'] = data['date'].apply(lambda dt: int(dt.weekday() < 5))
data['year'] = data['date'].apply(lambda dt: dt.year % 2013)
data['month'] = data['date'].apply(lambda dt: dt.month)
data.head()
И наконец-таки сформируем выборку и ответы
X, y = data.drop(['target', 'date'], axis=1), data['target']
full_df = pd.concat([X, y], axis=1)
full_df.rename(mapper={0: 'target'}, axis=1, inplace=True)
Теперь, когда мы сформировали обучабщую выборку, давайте посмотрим на матрицу корреляций.
plt.subplots(figsize=(14,10))
sns.heatmap(full_df.corr(), cmap="BuPu");
Как и следовало ожидать, у погодных признаков одной категории высокая корреляция (напрмер, максимальная, средняя и минимальные температуры будут коррелировать). Это может спровоцировать линейные моделей на огромные весовые коэффициенты, но мы всегда можем избавиться от лишних признаков или воспользоваться регуляризацией.
Так же можно заметить, что целевой коррелирует с признаками "Будни" и "Выходные", причем у признака "Будни" корреляция выше. Ещё имеем отрицательную корреляцию с признаком "Rain-Thunderstorm" (т.е. дождь с грозой), что тоже вполне логично.
Соберём все данные и отложим последние три месяца для теста качества модели.
train_size = X.shape[0] - 90
X_train, y_train = X[:train_size], y[:train_size]
X_test, y_test = X[train_size:], y[train_size:]
Как всем известно, типичные метрики в задачах регрессии это средняя абсолютная (MAE) и среднеквадратичная (MSE) ошибки. MSE сильнее штрафует за большие отклонения и хорошо подходит для контроля качества во время обучения, но не позволяет сделать выводы о качестве полученного решения.
Ещё есть R2
мера:
$$\large
R^2 = 1 - \frac{\sum_{i=1}^l (a(x_i) - y_i^2)}{\sum_{i=1}^l (y_i - \hat{y})^2}$$
Чем ближе R2
к единице, тем лучше модель объясняет данные.
Если R2
близка к нулю, предсказания близки к константным.
Фактически данная мера - нормированная среднеквадратичная ошибка. Можно сказать получаем интерпретируемую MSE. Её мы и будем использовать в данной задаче.
def score_model(clf, title):
tscv = TimeSeriesSplit(n_splits=10)
cv_res = np.median(cross_val_score(clf, X_train, y_train, cv=tscv, scoring='r2'))
print(title, "R2: {0:.3f}".format(cv_res))
Попробуем обучить несколько линейных и "деревянных" моделей со стандартными параметрами(пусть это и не совсем честно), проведём кросс-валидацию по TimeSeriesSplit разбиению и проверим качество на отложенной выборке для того, что бы определиться с моделью.
Использование TimeSeriesSplit позволяет получить последовательные обучающие выборки из прошлого и тестовую из будущего, что приближает задачу к реальной. Первое разбиение TimeSeriesSplit будет очень маленьким и качество обучения на нём, соответсвенно, будет негативно сказывается на среднем качестве по кросс-валидации, поэтому вместо среднего будем рассматривать медиану.
models = [
(LinearRegression(), "Линейная регрессия"),
(Ridge(random_state=17), "Регрессия с регуляризатором L2"),
(Lasso(random_state=17), "Регрессия с регуляризатором L1"),
(RandomForestRegressor(random_state=17), "Случайный лес"),
(GradientBoostingRegressor(random_state=17), "Градиентный бустинг(sklearn)"),
(XGBRegressor(random_state=17), "Градиентный бустинг(xgb)")
]
for pair in models:
score_model(*pair)
Модель линейной регрессии показала самый слабый результат. Это можно объяснить как следствие мультиколлинеарности. Регрессия с регуляризацией, показала схожий результат. Lasso и Ridge могут лучше, если подобрать параметр регуляризации. Случайный лес и бустинги имеют сравнимые результаты.
Выберем в качестве основной модели xgb (у него всё таки лучший скор на кросс валидации) и подберём параметры.
xgb = XGBRegressor(random_state=17)
tscv = TimeSeriesSplit(n_splits=10)
parameters = {'n_estimators': range(10, 311, 30),
'max_depth': [3, 9, 15, 45],
}
grid = GridSearchCV(xgb, parameters, n_jobs=5,
cv=tscv, scoring='r2', refit=True)
grid.fit(X_train, y_train)
grid.best_params_
xgb = XGBRegressor(random_state=17, **grid.best_params_)
param_range = list(range(10, 200, 10))
tscv = TimeSeriesSplit(n_splits=10)
train_scores, test_scores = validation_curve(
xgb, X_train, y_train,
param_name="n_estimators", param_range=param_range,
cv=tscv, scoring='r2', n_jobs=-1)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)
test_scores_std = np.std(test_scores, axis=1)
layout = go.Layout(title="Кривые валидации",
xaxis=dict(title="n_estimators"),# type='log'),
yaxis=dict(title="Score"))
lw = 2
p1 = go.Scatter(x=param_range, y=train_scores_mean,
name="Training score",
mode='lines',
line=dict(color="orange", width=lw))
p2 = go.Scatter(x=param_range, y=train_scores_mean - train_scores_std,
mode='lines', showlegend=False,
line=dict(color="orange", width=1))
p3 = go.Scatter(x=param_range, y=train_scores_mean + train_scores_std,
mode='lines', showlegend=False,
line=dict(color="orange", width=1),
fill='tonexty')
p4 = go.Scatter(x=param_range, y=test_scores_mean,
name="Cross-validation score",
mode='lines',
line=dict(color="navy", width=lw))
p5 = go.Scatter(x=param_range, y=test_scores_mean - test_scores_std,
mode='lines', showlegend=False,
line=dict(color="navy", width=1))
p6 = go.Scatter(x=param_range, y=test_scores_mean + test_scores_std,
mode='lines', showlegend=False,
line=dict(color="navy", width=1),
fill='tonexty')
fig = go.Figure(data=[p2, p3, p5, p6, p1, p4], layout=layout)
iplot(fig, show_link=False)
После 50 деревьев увеличение сложности по данному параметру не улучшает качество модели, но кривые валидации всё ещё находятся достаточно далеко друг от друга. Следовательно, увеличение сложности модели, например, по другим параметрам, вероятно повысит качество модели.
def plot_learning_curve(estimator, X, y, cv=None,
n_jobs=1, train_sizes=np.linspace(.1, 1.0, 5)):
train_sizes, train_scores, test_scores = learning_curve(
estimator, X, y,
cv=cv, n_jobs=n_jobs, train_sizes=train_sizes,
scoring='r2', random_state=17
)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)
test_scores_std = np.std(test_scores, axis=1)
layout = go.Layout(title="Кривые обучения")
p1 = go.Scatter(x=train_sizes, y=test_scores_mean + test_scores_std,
mode='lines',
line=dict(color="green", width=1),
showlegend=False)
p2 = go.Scatter(x=train_sizes, y=test_scores_mean - test_scores_std,
mode='lines',
line=dict(color="green", width=1),
showlegend=False, fill='tonexty')
p3 = go.Scatter(x=train_sizes, y=train_scores_mean + train_scores_std,
mode='lines',
line=dict(color="red", width=1),
showlegend=False)
p4 = go.Scatter(x=train_sizes, y=train_scores_mean - train_scores_std,
mode='lines',
line=dict(color="red", width=1),
showlegend=False, fill='tonexty')
p5 = go.Scatter(x=train_sizes, y=train_scores_mean,
marker=dict(color='red'),
name="Training score", showlegend=True)
p6 = go.Scatter(x=train_sizes, y=test_scores_mean,
marker=dict(color='green'),
name="Cross-validation score", showlegend=True)
fig = go.Figure(data=[p1, p2, p3, p4, p5, p6], layout=layout)
iplot(fig, show_link=False)
xgb = XGBRegressor(random_state=17, **grid.best_params_)
tscv = TimeSeriesSplit(n_splits=10)
plot_learning_curve(xgb, X_train, y_train, cv=tscv)
Кривые обучения не сошлись, информации недостаточно и добавление новых данных может улучшить результат.
При формировании обучающей выборки, мы сделали отложенную выборку из трёх последних месяцев, проверим насколько хороша наша модель.
xgb = XGBRegressor(random_state=17, **grid.best_params_)
xgb.fit(X_train, y_train)
y_pred = xgb.predict(X_test)
print("R2: {0:.3f}, MAE: {1:.3f}".format(r2_score(y_test, y_pred),
mean_absolute_error(y_test, y_pred)))
Мы получили достаточно неплохие результаты R2 - 0.84, т.е. модель хорошо объясняет данные, MAE - 108 т.е. модель ошибается в среднем всего на 100 велосипедов, тоже приемлимо учитывая масштабы целевой переменной.
y_train.median(), y_test.median()
Так же можно посмотреть на важность признаков
importances = xgb.feature_importances_
indices = np.argsort(importances)
plt.subplots(figsize=(14,10))
plt.title('Важность признаков')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), X_train.columns[indices])
plt.show()
Получили интересную модель предсказывающую количество поездок на велосипедах по погодным условиям. Самые важные, по мнению модели, признаки это:
Качество модели на тестовой выборке по метрике R2 составляет 0.838, что означает хорошую обобщающую способность.
Возможные улучшения: