Физтех.Статистика
Информация не актуальна. Сайт переехал на miptstats.github.io
Скачать ipynb
Введение в анализ данных¶
Парадокс времени ожидания: а что в реальной жизни?¶
В предыдущем примере мы посмотрели на то, какие эффекты возникают в случае если автобусы прибывают на остановку с одинаковой интенсивностью независимо друг от друга. Выполняется ли это в реальности?
Исследуем данные о запланированном и фактическом времени прибытия для автобусов RapidRide маршрутов C, D и E на автобусной остановке 3rd&Pike в центре Сиэтла. Ниже приведена схема общественного транспорта Сиэтла, на которой красными линиями отмечены рассматриваемые автобусные маршруты.
Сами данные и идея исследования взяты из статьи на Хабре.
import numpy as np
import pandas as pd
from random import choices, shuffle
import scipy.stats as sps
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='whitegrid', font_scale=1.3, palette='Set2')
Загрузим данные с помощью библиотеки pandas, с которой вы уже должны были познакомиться. Напечатаем также несколько первых строк таблицы. В данных есть следующая информация:
OPD_DATE— дата события;VEHICLE_ID— номер транспортного средства;RTE— номер маршрута;DIR— направление движения (на север или на юг);TRIP_ID— идентификатор рейса;STOP_ID— идентификатор остановки;STOP_NAME— наименование остановки;SCH_STOP_TM— время прибытия по расписанию (schedule);ACT_STOP_TM— время прибытия по факту (actual);
df = pd.read_csv('arrival_times.csv')
df.head()
1. Обработка и первичный анализ данных¶
В данных всего 39 157 записей
df.shape
Проверим, есть ли пропуски
df.isna()
На первый взгляд не видно. Посчитаем долю пропусков по столбцам, усреднив логические константы. В данном случае сначала произошло преобразование к 0 и 1.
df.isna().mean()
Теперь уже видно, что для 0.6% рейсов не задано время прибытия.
Посмотрим, что это за рейсы. Для этого сначала с помощью функции any определим для каждого рейса, есть для него хотя бы один пропуск.
df.isna().any(axis=1)
Теперь мы можем извлеть все строки таблицы с пропусками.
df[df.isna().any(axis=1)]
Давайте удали пропуски из данных, чтобы они не мешали нашему дальнейшему анализу. Параметры функции предписывают удалять строки, в которых есть хотя бы один пропуск в любой из колонок.
df = df.dropna(axis=0, how='any')
Посмотрим на диаппазон дат в датасете. Всего у нас есть данные за 2 месяца — с 26 марта 2016 по 27 мая 2016.
df['OPD_DATE'].min(), df['OPD_DATE'].max()
Уникальные идентификаторы автобусных маршрутов
df['RTE'].unique()
Направления движения
df['DIR'].unique()
Идентификаторы автобусных остановок
df['STOP_ID'].unique()
Как видим остановки две. Если мы составим с помощью функции `pd.crosstab` таблицу, в которой показано, сколько рейсов было на каждой остановке по каждому из направлений, то мы сможем догадаться, что два идентификатора означаю одну и ту же остановку по разные стороны движения.
pd.crosstab(columns=df['DIR'], index=df['STOP_ID'])
Также с помощью этой функции мы можем составить таблицу с количеством автобусов каждого маршрута за каждые сутки
daily_counts = pd.crosstab(columns=df['RTE'], index=df['OPD_DATE'])
daily_counts
И построить с ее помощью график автобусов за день по маршрутам. На этом графике мы четко наблюдаем некоторую периодичность, которую называют недельной сезонностью. Действительно, мы видим, что 7-дневные интервалы, из которых в течении пяти дней наблюдается достаточно большое количество автобусов по каждому из маршрутов, а в остальные два дня — существенно меньше.
В целом стоит отметить, что исследование сезонностей в данных — важная часть работы с временными данными.
with sns.axes_style('darkgrid'):
daily_counts.plot(figsize=(15, 5), lw=2);
plt.ylabel('Количество автобусов')
Для удобства дальнейшего анализа объединим дату и время в одну метку времени.
df['scheduled'] = pd.to_datetime(df['OPD_DATE'] + ' ' + df['SCH_STOP_TM'])
df['actual'] = pd.to_datetime(df['OPD_DATE'] + ' ' + df['ACT_STOP_TM'])
Однако, такой операцией мы могли привнести ошибки в данные.
Давайте посчитаем для каждого автобуса время между прибытием по расписанию и фактическим временем прибытия.
minute = np.timedelta64(1, 'm')
hour = 60 * minute
diff_hrs = (df['actual'] - df['scheduled']) / hour
Посмотрим на гистограмму
plt.figure(figsize=(12, 5))
plt.hist(diff_hrs, bins=48, log=True)
plt.xticks(range(-24, 25, 3))
plt.xlabel('Отклонения от расп., час')
plt.ylabel('Кол-во автобусов');
Оказывается, что несколько автобусов приезжали на сутки раньше, а еще несколько — на сутки позже. Как такое могло произойти?
Вспомним, что мы искусственно склеили дату с временем суток. У нас имеется две отметки времени суток для каждого рейса, но дата одна, и она соответствует времени прибытия по расписанию. Таким образом, возможны две следующие ситуации.
- Время прибытия по расписанию — сразу после полуночи, а автобус пришел немного раньше. Тем самым при склейке даты и времени мы получили, что как будто он пришел на сутки позже.
- Время прибытия по расписанию — непосредственно до полуночи, а автобус немного опоздал. И при склейке даты и времени мы получили, что как будто он пришел на сутки раньше.
Исправим фактическое время прибытия руками для тех рейсов, где отклонение по модулю слишком большое.
df.loc[diff_hrs > 20, 'actual'] -= 24 * hour
df.loc[diff_hrs < -20, 'actual'] += 24 * hour
И теперь посчитаем время отклонения факта от расписания в минутах
df['minutes_late'] = (df['actual'] - df['scheduled']) / minute
Переобозначим технические идентификаторы маршрутов на те, которые видят пассажиры. Также переобозначим направления движения.
df['route'] = df['RTE'].replace({673: 'C', 674: 'D', 675: 'E'}).astype('category')
df['direction'] = df['DIR'].replace({'N': 'север', 'S': 'юг'}).astype('category')
Оставим только интересующие нас столбцы
df = df[['route', 'direction', 'scheduled', 'actual', 'minutes_late']]
И посмотрим на данные
df.head()
2. Исследование отклонения от расписания¶
Для каждого направления и каждого маршрута визуализируем гистограмму отклонения от расписания в минутах.
Визуализацию выполним в виде таблицы с помощью класса `FacetGrid`. Работа функции несколько похожа на функцию `pd.pivot_table`, только вместо применения аггрегирующией функции рисуется график. Функция построения каждого графика передается в метод map.
g = sns.FacetGrid(df, row="direction", col="route", height=3, aspect=2)
g.map(plt.hist, "minutes_late", bins=np.arange(-10, 15))
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('Отклонения от расп., мин.', 'Кол-во автобусов');
Заметим, что автобусы ближе к графику в начале маршрута и больше отклоняются от него к концу. Если внимательно посмотреть на схему маршрутов, то можно увидеть, что исследуемая остановка находится на северной части маршрута С, и на южных частях машрутов D и Е. Тем самым, при следовании по маршруту C на юг наша остановка близка к началу маршрута, тем самым отклонения от расписания в основном небольшие. Аналогично с маршрутами D и E при движении на север. В остальных случаях наблюдается больший разброс отклонений.
Однако нам хотелось бы еще посмотреть на интервалы между автобусами. Сначала посчитаем интервалы между автобусами по расписанию и фактические.
def compute_headway(times):
""" Вычисляет интервалы между соседними моментами времени в минутах. """
minute = np.timedelta64(1, 'm')
return times.sort_values().diff() / minute
# сгруппируем по маршруту и по направлению
grouped = df.groupby(['route', 'direction'])
# для каждой группы посчитаем интервалы определенной выше функцией
df['actual_interval'] = grouped['actual'].transform(compute_headway)
df['scheduled_interval'] = grouped['scheduled'].transform(compute_headway)
Визуализируем фактические интервалы между автобусами
g = sns.FacetGrid(df.dropna(), row="direction", col="route", height=3, aspect=2)
g.map(plt.hist, "actual_interval", bins=np.arange(40))
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('Актуальные интервалы, мин.', 'Кол-во автобусов');
Видим, что фактические интервалы даже по каждому отдельному маршруту, видимо, не подчиняются экспоненциальному распределению.
Однако, не все так просто. В прошлый раз мы рассматривали случай, при котором автобусы приходят на остановку независимо друг от друга с одинаковой интенсивностью. Но это не так, для этого достаточно посмотреть на графики временных рядов выше. Как мы ранее заметили, в выходные дни автобусы ходят реже. Кроме того, наверняка в час пик автобусы ходят чаще, чем в вечерне-ночное время.
Посмотрим на интервалы между автобусами по расписанию
g = sns.FacetGrid(df.dropna(), row="direction", col="route", height=3, aspect=2)
g.map(plt.hist, "scheduled_interval", bins=np.arange(40))
g.set_titles('{col_name} {row_name}')
g.set_axis_labels('Интервалы по расписанию, мин.', 'Кол-во автобусов');
Тут мы уже видим, что по расписанию есть много разных значений запланированных интервалов между автобусами.
3. Исследование для однородных ожидаемых интервалов¶
Тем не менее, есть частые интервалы по расписанию: 10, 12 и 15 минут. Например, есть почти 2000 автобусов в северную сторону с запланированным интервалом в 10 минут. Давайте посмотрим внимательнее на эти интервалы.
Составим условие выбора строк из таблицы
df['scheduled_interval'].isin([10, 12, 15])
Составим подтаблицу с помощью данного условия
subset = df[df['scheduled_interval'].isin([10, 12, 15])]
Сгруппируем по маршруту, по направлению и по интервалам среди выбранных интервалов.
grouped = subset.groupby(['route', 'direction', 'scheduled_interval'])
def stack_sequence(data):
""" Перераспределение времени так, как будто события происходят последовательно. """
# сортируем по ожидаемому времени прибытия
data = data.sort_values('scheduled')
# переопределим время прибытия так, как будто они шли подряд
data['scheduled'] = data['scheduled_interval'].cumsum()
# соответствующе переопределим фактическое время на основе имеющегося отклонения
data['actual'] = data['scheduled'] + data['minutes_late']
# посчитаем фактические интервалы по скорректированному фактическому времени
data['actual_interval'] = data['actual'].sort_values().diff()
return data
Применим эту функцию к нашим группам
sequenced = grouped.apply(stack_sequence).reset_index(drop=True)
sequenced.head()
Теперь визуализируем все отклонения от расписания в рамках выбранного интервала.
for route in ['C', 'D', 'E']:
g = sns.FacetGrid(sequenced.query(f"route == '{route}'"),
row="direction", col="scheduled_interval", height=3, aspect=2)
g.map(plt.hist, "actual_interval", bins=np.arange(40) + 0.5)
g.set_titles('{row_name} ({col_name:.0f} мин.)')
g.set_axis_labels('Актуальные интервалы, мин.', 'Кол-во автобусов');
g.fig.suptitle(f'{route} line', y=1.05, fontsize=28)
Видим, что для каждого маршрута и направления при фиксированном интервале распределение фактических интервалов близко к нормальному. Оно достигает максимума около запланированного интервала и имеет стандартное отклонение (корень из дисперсии), которое меньше в начале маршрута (на юг для C, на север для D/E) и больше в конце.
Также видно, что фактические интервалы прибытия определенно не соответствуют экспоненциальному распределению, что является основным предположением, на котором основан парадокс времени ожидания.
4. Симуляция пассажиров¶
Теперь оценим, сколько в среднем придется ждать пассажиру автобуса при данных интервалах движения. Для этого оформим код с предыдущего ноутбука в специальную функцию.
def simulate_wait_times(bus_arrival_times, n_passengers=1000000):
"""
Моделирование времени ожидание пассажиров.
bus_arrival_times -- моменты времени прибытия автобусов.
n_passengers -- количество пассажиров для семплирования.
"""
bus_arrival_times = np.array(bus_arrival_times)
# сгенерируем для каждого пассажира время его прибытия на остановку
passenger_times = sps.uniform(scale=bus_arrival_times.max()).rvs(size=n_passengers)
# найдем время прибытия следующего автобуса поиском по отсортированному массиву
i = np.searchsorted(bus_arrival_times, passenger_times, side='right')
# вычислим интервал ожидания
wait_times = bus_arrival_times[i] - passenger_times
return wait_times
Произведем семплирование и посчитаем интервалы вместе с наиболее вероятным диаппазоном, считая, что мы работаем с нормальным распределением. Как мы помним, интервал, определяемый среднем и двумя стандартными отклонениями содержит чуть более 95% вероятностной массы.
# сгруппируем по маршрутам, направлениям и интервалам
grouped = sequenced.groupby(['route', 'direction', 'scheduled_interval'])
# применим семплирование
simulated = grouped['actual'].apply(simulate_wait_times)
# посчитаем и напечатаем интервалы
simulated = simulated.apply(lambda times: "{0:.1f} +/- {1:.1f}".format(times.mean(), 2*times.std()))
pd.DataFrame(simulated)
Видим, что среднее время ожидания, возможно, на минуту или две больше половины интервала расписанию, но не равно этому интервалу, как подразумевает парадокс времени ожидания. Другими словами, парадокс инспекции подтвержден, но парадокс времени ожидания не соответствует действительности.
Вывод.
Мы подтвердили, что в реальном мире автобусные маршруты подчиняются некоторой разновидности парадокса инспекции, однако основное предположение, лежащее в основе парадокса времени ожидания, — что прибытие автобусов происходит независимо друг от друга и с одинаковой интенсивностью — не является действительным.
На самом же деле в хорошо управляемой системе общественного транспорта есть специально структурированные расписания, чтобы избежать такого поведения: автобусы не начинают свои маршруты в случайное время в течение дня, а стартуют по расписанию, выбранному для наиболее эффективной перевозки пассажиров.