Введение в анализ данных

Время ожидания автобуса

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")

1. Когда придет мой автобус? Или каково среднее время ожидания автобуса.

Вы приходите на автобусную остановку. Согласно расписанию автобус ходит каждые 10 минут. Вы засекаете время, и получается, что автобусы обычно приходят через 9-11 минут.

Почему так не везет?

Казалось бы, если автобус ходит каждые 10 минут, то в среднем его нужно ждать 5 минут. Однако, не все так просто. В действительности автобусы приезжают не точно по расписанию, а случайным образом. Оказывается, справедливо следующее утверждение.

Если автобусы приходят с одинаковой интенсивностью в 10 минут и независимо друг от друга, то среднее время ожидания составляет 10 минут.

Это называют парадоксом времени ожидания. Далее мы попробуем экспериментально проверить это утверждение.

Для точности вычислений возьмем достаточно большую выборку — миллион автобусов, которые приходят с интенсивностью раз в 10 минут и сгенерируем интервалы прибытия между автобусами.

count = 1000000  # количество автобусов
tau = 10  # средний интервал между автобусами

bus_arrival_times = sps.uniform(scale=count * tau).rvs(size=count)
bus_arrival_times = np.sort(bus_arrival_times)

Посмотрим на гистограмму интервалов прибытия автобусов, сравнивая ее с графиком плотность равномерного распределения.

grid = np.linspace(-10, count * tau + 10, 1010)

plt.figure(figsize=(8, 5))
plt.hist(bus_arrival_times, bins=50, density=True, label="Данные")
plt.plot(
    grid,
    sps.uniform(scale=count * tau).pdf(grid),
    lw=5,
    alpha=0.5,
    label="Плотность $U[0, max]$",
    color="#FF3300",
)
plt.xlabel("Время прибытия автобуса")
plt.legend()
plt.show()

Прежде чем перейти к ответу на наш исходный вопрос, посмотрим, сколько в среднем автобусов приезжает на остановку в течение часа? Поскольку у нас может не получится ровное количество часов, посчитаем для первых 100 000 часов.

hours_count = 100000  # количество часов

# количество автобусов за каждый интервал длиной в 1 час
hist = np.histogram(bus_arrival_times, bins=60 * np.arange(hours_count))[0]
# количество интервалов, за которые приехало одинаковое количество автобусов
x, y = np.unique(hist, return_counts=True)

Посмотрим на гистограмму распределения количества автобусов за час и сравним его с графиком дискретной плотности пуассоновского распределения с параметром 6. Число 6 есть теоретически ожидаемое количество автобусов в час если их интенсивность движения есть в среднем 1 автобус за 10 минут.

grid = np.arange(25)

plt.figure(figsize=(8, 5))
plt.bar(x, y / y.sum(), label="Данные")
plt.plot(
    grid,
    sps.poisson(mu=60 / tau).pmf(grid),
    marker="o",
    ms=10,
    lw=5,
    alpha=0.5,
    label="Дискр.\nплотность $Pois(60/\\tau)$",
    color="#FF3300",
)
plt.xlabel("Количество автобусов в час")
plt.legend()
plt.show()

Теперь рассмотрим средние длины интервалов между прибытиями автобусов. Посчитаем также средний интервал. Как видим, он составляет примерно 10 минут.

intervals = np.diff(bus_arrival_times)
intervals.mean()

np.float64(9.999991940498665)

Посмотрим на гистограмму длин интервалов между прибытиями автобусов и сравним ее с плотностью экспоненциального распределения.

grid = np.linspace(-1, 70, 1000)

plt.figure(figsize=(8, 5))
plt.hist(intervals, bins=50, density=True, label="Данные", range=(0, 70))
plt.plot(
    grid,
    sps.expon(scale=tau).pdf(grid),
    lw=5,
    alpha=0.5,
    label="Плотность $Exp(\\tau)$",
    color="#FF3300",
)
plt.xlabel("Время между прибытиями автобусов")
plt.legend()
plt.show()

Проведем опрос миллиона пассажиров с целью выяснить, сколько времени они ждали автобус.

n_passengers = 1000000  # количество пассажиров

# сгенерируем для каждого пассажира время его прибытия на остановку
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

wait_times.mean()

np.float64(10.010549715321439)

Посмотрим на гистограмму времени ожидания прибытия автобуса и сравним ее с плотностью ТОГО ЖЕ САМОГО экспоненциального распределения.

grid = np.linspace(-1, 70, 1000)

plt.figure(figsize=(8, 5))
plt.hist(wait_times, bins=50, density=True, label="Данные", range=(0, 70))
plt.plot(
    grid,
    sps.expon(scale=tau).pdf(grid),
    lw=5,
    alpha=0.5,
    label="Плотность $Exp(\\tau)$",
    color="#FF3300",
)
plt.xlabel("Время ожидания автобуса")
plt.legend()
plt.show()

Оказывается, оставшееся время ожидания автобуса не зависит от того времени, которое уже прошло с момента прибытия предыдущего автобуса. Это свойство называется свойством отсутствия памяти. Среди всех абсолютно непрерывных распределений таким свойством обладает только экспоненциальное распределение. Среди всех дискретных распредений — только геометрическое.

Если мы считаем, что у нас есть несколько маршрутов автобусов, причем все автобусы приходят независимо друг от друга с одинаковой интенсивностью, то количество автобусов, которое придется пропустить, имеет геометрическое распределение.

А сколько времени проходит между каждым пятым прибывающим автобусом? Посчитаем эти интервалы времени.

intervals = np.diff(bus_arrival_times[::5])
intervals.mean()

np.float64(50.00007057551471)

grid = np.linspace(-1, 180, 1000)

plt.figure(figsize=(8, 5))
plt.hist(intervals, bins=50, density=True, label="Данные", range=(0, 180))
plt.plot(
    grid,
    sps.gamma(a=5, scale=tau).pdf(grid),
    lw=5,
    alpha=0.5,
    label="Плотность $\\Gamma(5, \\tau)$",
    color="#FF3300",
)
plt.xlabel("Время между прибытиями каждого 5-го автобуса")
plt.legend()
plt.show()

Этот пример также показывает, что экпоненциальное распределение является частным случаем гамма-распределения.

2. Парадокс времени ожидания: а что в реальной жизни?

В предыдущем примере мы посмотрели на то, какие эффекты возникают в случае если автобусы прибывают на остановку с одинаковой интенсивностью независимо друг от друга. Выполняется ли это в реальности?

Исследуем данные о запланированном и фактическом времени прибытия для автобусов RapidRide маршрутов C, D и E на автобусной остановке 3rd&Pike в центре Сиэтла. Ниже приведена схема общественного транспорта Сиэтла, на которой красными линиями отмечены рассматриваемые автобусные маршруты.

Сами данные и идея исследования взяты из статьи на Хабре.

Загрузим данные с помощью библиотеки 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()

	OPD_DATE	VEHICLE_ID	RTE	DIR	TRIP_ID	STOP_ID	STOP_NAME	SCH_STOP_TM	ACT_STOP_TM
0	2016-03-26	6201	673	S	30908177	431	3RD AVE & PIKE ST (431)	01:11:57	01:13:19
1	2016-03-26	6201	673	S	30908033	431	3RD AVE & PIKE ST (431)	23:19:57	23:16:13
2	2016-03-26	6201	673	S	30908028	431	3RD AVE & PIKE ST (431)	21:19:57	21:18:46
3	2016-03-26	6201	673	S	30908019	431	3RD AVE & PIKE ST (431)	19:04:57	19:01:49
4	2016-03-26	6201	673	S	30908252	431	3RD AVE & PIKE ST (431)	16:42:57	16:42:39

Обработка и первичный анализ данных

В данных всего 39 157 записей.

df.shape

(39157, 9)

Проверим, есть ли пропуски.

df.isna()

	OPD_DATE	VEHICLE_ID	RTE	DIR	TRIP_ID	STOP_ID	STOP_NAME	SCH_STOP_TM	ACT_STOP_TM
0	False	False	False	False	False	False	False	False	False
1	False	False	False	False	False	False	False	False	False
2	False	False	False	False	False	False	False	False	False
3	False	False	False	False	False	False	False	False	False
4	False	False	False	False	False	False	False	False	False
...	...	...	...	...	...	...	...	...	...
39152	False	False	False	False	False	False	False	False	False
39153	False	False	False	False	False	False	False	False	False
39154	False	False	False	False	False	False	False	False	False
39155	False	False	False	False	False	False	False	False	False
39156	False	False	False	False	False	False	False	False	False

39157 rows × 9 columns

На первый взгляд не видно. Посчитаем долю пропусков по столбцам, усреднив логические константы. В данном случае сначала произошло преобразование к 0 и 1.

df.isna().mean()

OPD_DATE       0.000000
VEHICLE_ID     0.000000
RTE            0.000000
DIR            0.000000
TRIP_ID        0.000000
STOP_ID        0.000000
STOP_NAME      0.000000
SCH_STOP_TM    0.006129
ACT_STOP_TM    0.000000
dtype: float64

Теперь уже видно, что для 0.6% рейсов не задано время прибытия.

Посмотрим, что это за рейсы. Для этого сначала с помощью функции any определим для каждого рейса, есть для него хотя бы один пропуск.

df.isna().any(axis=1)

0        False
1        False
2        False
3        False
4        False
         ...  
39152    False
39153    False
39154    False
39155    False
39156    False
Length: 39157, dtype: bool

Теперь мы можем извлеть все строки таблицы с пропусками.

df[df.isna().any(axis=1)]

	OPD_DATE	VEHICLE_ID	RTE	DIR	TRIP_ID	STOP_ID	STOP_NAME	SCH_STOP_TM	ACT_STOP_TM
383	2016-04-01	6206	673	S	0	431	3RD AVE & PIKE ST (431)	NaN	16:00:57
471	2016-03-31	6207	673	S	0	431	3RD AVE & PIKE ST (431)	NaN	15:28:05
588	2016-03-30	6208	673	S	0	431	3RD AVE & PIKE ST (431)	NaN	14:52:36
599	2016-03-30	6200	673	S	0	431	3RD AVE & PIKE ST (431)	NaN	15:12:06
685	2016-03-29	6213	673	S	0	431	3RD AVE & PIKE ST (431)	NaN	17:56:29
...	...	...	...	...	...	...	...	...	...
37327	2016-05-10	6215	674	N	0	578	3RD AVE & PIKE ST (578)	NaN	15:25:47
37328	2016-05-10	6215	674	N	0	578	3RD AVE & PIKE ST (578)	NaN	17:44:05
37570	2016-05-12	6218	674	N	0	578	3RD AVE & PIKE ST (578)	NaN	18:27:55
37670	2016-05-13	6212	674	N	0	578	3RD AVE & PIKE ST (578)	NaN	19:04:25
39107	2016-05-27	6099	674	N	0	578	3RD AVE & PIKE ST (578)	NaN	14:51:34

240 rows × 9 columns

Давайте удалим пропуски из данных, чтобы они не мешали нашему дальнейшему анализу. Параметры функции предписывают удалять строки, в которых есть хотя бы один пропуск в любой из колонок.

df = df.dropna(axis=0, how="any")

Посмотрим на диапазон дат в датасете. Всего у нас есть данные за 2 месяца — с 26 марта 2016 по 27 мая 2016.

df["OPD_DATE"].min(), df["OPD_DATE"].max()

('2016-03-26', '2016-05-27')

Уникальные идентификаторы автобусных маршрутов.

df["RTE"].unique()

array([673, 675, 674])

Направления движения.

df["DIR"].unique()

array(['S', 'N'], dtype=object)

Идентификаторы автобусных остановок.

df["STOP_ID"].unique()

array([431, 578])

Как видим остановки две. Если мы составим с помощью функции pd.crosstab таблицу, в которой показано, сколько рейсов было на каждой остановке по каждому из направлений, то мы сможем догадаться, что два идентификатора означают одну и ту же остановку по разные стороны движения.

pd.crosstab(columns=df["DIR"], index=df["STOP_ID"])

DIR	N	S
STOP_ID
431	0	19318
578	19599	0

Также с помощью этой функции мы можем составить таблицу с количеством автобусов каждого маршрута за каждые сутки.

daily_counts = pd.crosstab(columns=df["RTE"], index=df["OPD_DATE"])
daily_counts

RTE	673	674	675
OPD_DATE
2016-03-26	176	178	159
2016-03-27	153	153	140
2016-03-28	236	228	221
2016-03-29	233	231	221
2016-03-30	230	230	221
...	...	...	...
2016-05-23	230	230	221
2016-05-24	236	229	220
2016-05-25	235	231	224
2016-05-26	231	203	222
2016-05-27	231	231	223

63 rows × 3 columns

И построить с ее помощью график автобусов за день по маршрутам. На этом графике мы четко наблюдаем некоторую периодичность, которую называют недельной сезонностью. Действительно, мы видим, что 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()

	route	direction	scheduled	actual	minutes_late
0	C	юг	2016-03-26 01:11:57	2016-03-26 01:13:19	1.366667
1	C	юг	2016-03-26 23:19:57	2016-03-26 23:16:13	-3.733333
2	C	юг	2016-03-26 21:19:57	2016-03-26 21:18:46	-1.183333
3	C	юг	2016-03-26 19:04:57	2016-03-26 19:01:49	-3.133333
4	C	юг	2016-03-26 16:42:57	2016-03-26 16:42:39	-0.300000

Исследование отклонения от расписания

Для каждого направления и каждого маршрута визуализируем гистограмму отклонения от расписания в минутах.

Визуализацию выполним в виде таблицы с помощью класса 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)

<ipython-input-17-5654ec461668>:8: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  grouped = df.groupby(["route", "direction"])

Визуализируем фактические интервалы между автобусами.

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("Интервалы по расписанию, мин.", "Кол-во автобусов");

Тут мы уже видим, что по расписанию есть много разных значений запланированных интервалов между автобусами.

Исследование для однородных ожидаемых интервалов

Тем не менее, есть частые интервалы по расписанию: 10, 12 и 15 минут. Например, есть почти 2000 автобусов в северную сторону с запланированным интервалом в 10 минут. Давайте посмотрим внимательнее на эти интервалы.

Составим условие выбора строк из таблицы.

df["scheduled_interval"].isin([10, 12, 15])

	scheduled_interval
0	False
1	True
2	True
3	True
4	True
...	...
39152	True
39153	False
39154	False
39155	False
39156	False

39157 rows × 1 columns

dtype: bool

Составим подтаблицу с помощью данного условия.

subset = df[df["scheduled_interval"].isin([10, 12, 15])]

Сгруппируем по маршруту, по направлению и по интервалам среди выбранных интервалов.

grouped = subset.groupby(["route", "direction", "scheduled_interval"])

<ipython-input-23-52883d3df64d>:1: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  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()

<ipython-input-25-911cc68e75f3>:1: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  sequenced = grouped.apply(stack_sequence).reset_index(drop=True)

	route	direction	scheduled	actual	minutes_late	actual_interval	scheduled_interval
0	C	север	10.0	12.133333	2.133333	NaN	10.0
1	C	север	20.0	22.633333	2.633333	10.500000	10.0
2	C	север	30.0	34.350000	4.350000	11.716667	10.0
3	C	север	40.0	37.516667	-2.483333	3.166667	10.0
4	C	север	50.0	48.733333	-1.266667	11.216667	10.0

Теперь визуализируем все отклонения от расписания в рамках выбранного интервала.

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) и больше в конце.

Также видно, что фактические интервалы прибытия определенно не соответствуют экспоненциальному распределению, что является основным предположением, на котором основан парадокс времени ожидания.

Симуляция пассажиров

Теперь оценим, сколько в среднем придется ждать пассажиру автобуса при данных интервалах движения. Для этого оформим код из предыдущего примера в специальную функцию.

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

Произведем семплирование.

# сгруппируем по маршрутам, направлениям и интервалам
grouped = sequenced.groupby(["route", "direction", "scheduled_interval"])
# применим семплирование
simulated = grouped["actual"].apply(simulate_wait_times)

# посчитаем средние значения при сэмплировании
simulated = simulated.apply(
    lambda times: "{0:.1f}".format(times.mean())
)
pd.DataFrame(simulated)

<ipython-input-35-5999c72c3141>:2: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
  grouped = sequenced.groupby(["route", "direction", "scheduled_interval"])

			actual
route	direction	scheduled_interval
C	север	10.0	7.0
		12.0	7.3
		15.0	8.3
	юг	10.0	5.9
		12.0	6.7
		15.0	7.8
D	север	10.0	5.4
		12.0	6.2
		15.0	7.8
	юг	10.0	6.8
		12.0	7.4
		15.0	8.5
E	север	10.0	5.5
		12.0	6.2
		15.0	7.8
	юг	10.0	7.0
		12.0	7.9
		15.0	8.7

Видим, что среднее время ожидания, возможно, на минуту или две больше половины интервала расписанию, но не равно этому интервалу, как подразумевает парадокс времени ожидания. Другими словами, парадокс инспекции подтвержден, но парадокс времени ожидания не соответствует действительности.

Вывод.

Мы подтвердили, что в реальном мире автобусные маршруты подчиняются некоторой разновидности парадокса инспекции, однако основное предположение, лежащее в основе парадокса времени ожидания, — что прибытие автобусов происходит независимо друг от друга и с одинаковой интенсивностью — не является действительным.

На самом же деле в хорошо управляемой системе общественного транспорта есть специально структурированные расписания, чтобы избежать такого поведения: автобусы не начинают свои маршруты в случайное время в течение дня, а стартуют по расписанию, выбранному для наиболее эффективной перевозки пассажиров.