Вступление

Сегодня мы рассмотрим интересный пример - прогноз спроса на такси в Нью-Йорке. Проще говоря, водитель ожидает пассажиров в любом уголке города, и мы можем предсказать количество пикапов в этих регионах на следующие 10 минут. Так водитель будет знать, сколько пикапов он может ожидать за 10 минут. В этом случае я уменьшил MAPE (средняя абсолютная ошибка в процентах) до 0,02. Я использовал тройное экспоненциальное сглаживание, что очень эффективно при работе с данными временных рядов.

Для этого тематического исследования мы заимствовали данные из Комиссии по такси и лимузинам Нью-Йорка. Актуальные данные можно найти здесь. Компания NYC Taxi and Limousine собрала очень много данных обо всех поездках на такси в городе.

Вы можете найти исходный код этого кейса здесь.

1. Информация о такси:

Желтое такси: такси с желтым медальоном

Это знаменитые желтые такси Нью-Йорка, которые передвигаются исключительно по улицам. Количество такси ограничено конечным количеством медальонов, выдаваемых TLC. Чтобы получить доступ к этому способу передвижения, нужно стоять на улице и рукой ловить доступное такси. Пикапы заранее не оговорены.

Для автомобилей напрокат (FHV)

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

Зеленое такси: ливрея Street Hail (SHL)

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

Кредиты: Quora

В этом тематическом исследовании я рассматриваю только желтые такси для периода с января по март 2015 года и с января по март 2016 года.

2. Обзор проблем бизнеса / реального мира

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

3. Сопоставление с проблемой машинного обучения:

Очевидно, что это проблема регрессии для данных временных рядов. Перед нами стоят следующие задачи:

  1. Как разбить Нью-Йорк на регионы.
  2. Каждый регион Нью-Йорка должен быть разбит на 10-минутный интервал.

4. Показатели производительности.

  1. Средняя абсолютная ошибка в процентах.
  2. Среднеквадратичная ошибка.

5. Информация о данных

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

Мы можем видеть, что на Рис. 01 есть 19 различных функций, таких как долгота возврата, широта высадки, дата и время посадки, долгота широты посадки и так далее. Обратите внимание, что каждая строка в наборе данных - это поездка.

5.1 Очистка данных

Очистка данных - самый важный шаг в машинном обучении. Очевидно, что данные могут иметь выбросы и ошибки. Проведем односторонний анализ.

5.1.1. Широта и долгота посадки

В этом примере мы будем рассматривать только Нью-Йорк, Нью-Йорк ограничен координатами местоположения (широта, долгота) - (40,5774, -74,15) и (40,9176, -73,7004), поэтому любые координаты вне этих координат просто отбрасываются. Нас интересуют только пикапы из Нью-Йорка. Исходный код находится здесь.

Мы можем заметить, что есть некоторые точки, которые указывают на место посадки такси в Океане, что, очевидно, является необычной точкой, поэтому я уберу такие точки. Буду рассматривать только координаты (широта, долгота) - (40.5774, -74.15), (40.9176, -73.7004).

5.1.2. Широта и долгота окончания

Подобно pickup lat, long нам нужно отфильтровать точки данных, которые не находятся в пределах нашей определенной границы, если мы их обнаружим, я просто отброшу их. Исходный код находится здесь.

5.1.3. Продолжительность поездки

В нашем наборе данных некоторые точки оказались отрицательными, очевидная продолжительность поездки должна быть только положительной. В то время как некоторые точки показывают 50к минут поездки. Используя график box cox, мы можем визуализировать выбросы для функции продолжительности поездки и удалить их. Исходный код находится здесь.

5.1.4. Скорость

Точно так же мы использовали коробчатый график Кокса для визуализации данных. Некоторая точка имеет скорость 192857142 миль / час, что, очевидно, является ошибкой. Я сниму такие точки. Средняя скорость в Нью-Йорке составляет 12,45 миль / час. Исходный код находится здесь.

5.1.5. Расстояние поездки

Точно так же максимальная дальность поездки, показанная с использованием процентиля, составляет 258,9 миль. Можно сказать, что это выброс. Так что сниму такие точки. Исходный код находится здесь.

5.1.6. Общая стоимость проезда

Опять же, здесь также максимальная общая стоимость проезда показывает 3950612 долларов, что практически невозможно, поэтому я уберу эти баллы.

После удаления всех выбросов и ошибочных точек размер данных отфильтрован следующим образом:

Removing outliers in the month of Jan-2015
----
Number of pickup records =  12748986
Number of outlier coordinates lying outside NY boundaries: 293919
Number of outliers from trip times analysis: 23889
Number of outliers from trip distance analysis: 92597
Number of outliers from speed analysis: 24473
Number of outliers from fare analysis: 5275
Total outliers removed 377910
---
fraction of data points that remain after removing outliers in % 97.03576425607495

6. Подготовка данных

После удаления выбросов и ошибочных точек наша следующая задача - разделить весь Нью-Йорк на сегменты, а также разбить весь временной интервал на 10-минутные интервалы.

6.1 Кластеризация / сегментация

Здесь я буду использовать алгоритм K-Means, чтобы разбить город на регионы. Алгоритм Kmeans имеет свойство создавать кластеры одинакового размера. При сегментации мы должны позаботиться о том, чтобы размер кластера не был слишком маленьким или слишком большим. Итак, нам помогут два условия:

  1. Максимальное расстояние между кластерами должно быть менее 2 миль.
  2. Минимальное расстояние между кластерами между двумя кластерами должно быть ›0,5 мили.

Другое дело, что нам нужно беспокоиться о том, что нам нужно найти оптимальное число K в алгоритме Kmean. Очевидно, мы должны настроить гиперпараметры и найти оптимальное минимальное значение K. Исходный код находится здесь.

#trying different cluster sizes to choose the right K in K-means
coords = frame_with_durations_outliers_removed[['pickup_latitude', 'pickup_longitude']].values
neighbours=[]

def find_min_distance(cluster_centers, cluster_len):
    nice_points = 0
    wrong_points = 0
    less2 = []
    more2 = []
    min_dist=1000
    for i in range(0, cluster_len):
        nice_points = 0
        wrong_points = 0
        for j in range(0, cluster_len):
            if j!=i:
                distance = gpxpy.geo.haversine_distance(cluster_centers[i][0], cluster_centers[i][1],cluster_centers[j][0], cluster_centers[j][1])
                min_dist = min(min_dist,distance/(1.60934*1000))
                if (distance/(1.60934*1000)) <= 2:
                    nice_points +=1
                else:
                    wrong_points += 1
        less2.append(nice_points)
        more2.append(wrong_points)
    neighbours.append(less2)
    print ("On choosing a cluster size of ",cluster_len,"\nAvg. Number of Clusters within the vicinity (i.e. intercluster-distance < 2):",
           np.ceil(sum(less2)/len(less2)), "\nAvg. Number of Clusters outside the vicinity (i.e. intercluster-distance > 2):", np.ceil(sum(more2)/len(more2)),
           "\nMin inter-cluster distance = ",min_dist,"\n---")

def find_clusters(increment):
    kmeans = MiniBatchKMeans(n_clusters=increment, batch_size=10000,random_state=42).fit(coords)
    frame_with_durations_outliers_removed['pickup_cluster'] = kmeans.predict(frame_with_durations_outliers_removed[['pickup_latitude', 'pickup_longitude']])
    cluster_centers = kmeans.cluster_centers_
    cluster_len = len(cluster_centers)
    return cluster_centers, cluster_len

# we need to choose number of clusters so that, there are more number of cluster regions 
#that are close to any cluster center
# and make sure that the minimum inter cluster should not be very less
for increment in range(10, 100, 10):
    cluster_centers, cluster_len = find_clusters(increment)
    find_min_distance(cluster_centers, cluster_len)

Приведенный выше код находит гиперпараметр от (10 до 100), и мы после запуска видим, что K = 40 работает намного лучше. Итак, мы выберем K = 40. Итак, простыми словами, мы можем разделить Нью-Йорк на 40 регионов. Итак, у нас может быть 40 центроидов. Нанесение этих центроидов на график Нью-Йорка.

Мы можем видеть, что большинство центроидов расположено в районе Манхэттена, где густонаселенные районы Нью-Йорка и других регионов редки.

6.2 Временной интервал

У нас есть сегментированные регионы, теперь нам нужно сегментировать время по ячейкам времени. Как мы знаем, мы разделим время на 10-минутный интервал. 10 минут = 600 секунд. Мы конвертируем весь формат DateTime в формат временной метки Unix. Затем разделите его на 600. Итак, если вы считаете, что январь 2016 года, есть 24 часа, 31 день и 60 секунд (24 * 31 * 60/10 = 4464 бина). Аналогично для февраля 2015 года = 4176бинов (24 * 29 * 60/10).

Итак, теперь у нас есть две части, 40 кластеров и 4464 корзины по данным за январь 2016 года. На февраль 40 кластеров и 4176 бункеров. Для марта 40 кластеров и 4464 ячейки. Исходный код находится здесь.

7. Сглаживание

Есть некоторые бункеры, в которых нет пикапов. В полночь почти все бункеры не подбираются. Необязательно, но большая часть пикапа предсказывает значение 0. Этот 0 будет головной болью, когда мы будем выполнять расчетную часть, он может привести к ошибке деления на ноль. Вместо этого мы можем сопоставить все нулевые значения с соседними значениями. Предположим, что t, t1, t2 - интервалы времени. при t = 50 срабатываний, при t1 = 0 срабатываний и при t2 = 100 срабатываний, поэтому мы будем сглаживать как (100 + 50 + 0) / 3 = 50, поэтому t = 50, t1 = 50 и t2 = 50. Исходный код "здесь".

Теперь у нас есть 4464 тайм-бина с номерами пикапов в каждых 40 регионах. помещаем все тайм-бины в очередь. Затем мы должны обрабатывать нижеследующие граничные случаи,

Case 1: When we have the last/last few values are found to be missing,hence we have no right-limit here
Case 2: When we have the missing values between two known values
Case 3: When we have the first/first few values are found to be missing,hence we have no left-limit here.

Случай 1: Рассмотрим 4 временных интервала, t1, t2, t3, t4, где значение срабатывания равно 10,2,0,0. Отсутствуют правые значения. Т.е. в t1, t2 имеется 10,2 срабатывания, а в t3, t4 - нулевое срабатывание. поэтому мы возьмем среднее (10 + 2 + 0 + 0) / 4 = 3 и после сглаживания наводок (3,3,3,3).

Случай 2: аналогичные интервалы времени t1, t2, t3, t4, где значение срабатывания равно 10,0,0,2. Средние значения отсутствуют. т.е. t1 = 10, t2 = 0,3 = 0 и t4 = 2 срабатывания. аналогично после сглаживания данных мы можем получить (3,3,3,3).

Случай 3: тот же интервал времени t1, t2, t3, t4 и 0,0,2,10. Значения слева отсутствуют. Точно так же после сглаживания данных мы можем получить (3,3,3,3).

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

8. Временные ряды и преобразования Фурье

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

Если вы ничего не знаете о преобразовании Фурье, отметьте здесь. Но в некоторых словах можно использовать преобразование Фурье для разложения формы волны / сигнала. С этого момента я буду использовать волны, форму волны и сигнальное слово как синонимы. Посмотрите на приведенную выше форму волны повторяющейся природы. На самом деле она состоит из нескольких волн с разной частотой и амплитудой. Предположим, что синусоидальная волна s1 имеет некоторую амплитуду и частоту, синусоидальная волна s2 имеет некоторую амплитуду и частоту, а волна s3 также имеет некоторую амплитуду и частоту, тогда она объединит все s1, s2 и s3 и будет составлять синусоидальную волну или сложение всех трех волн. Назовем это "S".

Когда мы применяем преобразование Фурье поверх S-волны, оно разлагается, и мы можем получить обратно s1, s2 и s3. Результирующие s1, s2 и s3 будут в «частотной области». В задачах временных рядов преобразование Фурье может возвращать наиболее ценные характеристики. Помните, что в «частотной области» мы можем очень эффективно анализировать сигнал. Мы применили преобразование Фурье к нашему сигналу.

Мы можем наблюдать на изображении Fig02, что у нас есть данные временного ряда, которые имеют повторяющийся характер. В одном дне у нас 1440 минут (24 часа * 60 минут = 1440 минут). Мы делим его на 10-минутный интервал, так что волна fig02 повторяется через 144 единицы (1440 минут / 10 минут = 144 минуты). На Фиг03 показано преобразование Фурье 24-часовых данных.

На Рис03. первый самый высокий пик - это постоянная составляющая преобразования Фурье. Мы проигнорируем компонент постоянного тока. Второй пик - 1/144, что соответствует 24 часам, третий пик - 1/72, что соответствует 12 часам дня. Четвертый пик - 1/36, что соответствует 6 часам дня. скоро.

9. Базовая модель

У нас есть две стратегии для базовой модели.

  1. Соотношение Особенности
  2. Функция предыдущего значения

Предположим, что это 15 января 2016 года, 17.40, просто указано, что это текущая временная метка, и я хочу спрогнозировать спрос на такси на следующие 10 минут. Затем 15 января 2016, 17:40, мы будем называть его P_t 2016.

В то же время P_t 2015 - это не что иное, 15 января 2015, 17.40. с этого момента я буду использовать те же обозначения.

Также обратите внимание, что P_t или R_t представляют текущее значение.

P_t-1 или R_t-1 представляют прошлое значение.

P_t + 1 или R_t + 1 представляют будущее значение.

Как мы видим, для функций соотношения предположим, что текущая отметка времени - t. Итак, R_t = P_t 2016 / P_t 2015. Предположим, будущая временная метка t + 1, формулировка меняется следующим образом: R_t + 1 = P_t + 1 2016 / P_t +1 2015. где P_t + 1 2016 - (прогнозируемое значение) и P_t + 1 2015 г. - фактическая / доступная стоимость.

P_t + 1 2016 = R_t + 1 * P_t + 1 2015. где P_t + 1 2015 будет известен.

Предыдущая функция - это самый простой способ. Она будет учитывать прошлое значение метки времени. скажем, что текущая отметка времени - t, тогда будет использоваться t-1, t-2, t-3,…. данные.

9.1.1. Простая скользящая средняя

Это действительно просто формулировка, это простое среднее из предыдущих n значений Ratio features. Лучшее значение для n - 3. Я проделал некоторую настройку гиперпараметров для значения n и обнаружил, что 3 - лучшее значение для нашего случая. Исходный код находится здесь.

Точно так же для «функции предыдущего значения» мы взяли среднее значение из предыдущих n значений, после настройки параметра я обнаружил, что n = 1 работает лучше.

Примечание. В простой скользящей средней мы даем равный вес t-1, t-2, t-3, вместо этого можем ли мы попробовать другую схему весов?

9.1.2. Взвешенное скользящее среднее

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

9.1.3. Экспоненциально взвешенное скользящее среднее

Предположим, я хочу спрогнозировать значение R’_t при текущей метке времени и R_t-1 (прошлое фактическое значение) и R’_t-1 (прошлое прогнозируемое значение). Так что я могу использовать оба значения для прогнозов на будущее. В приведенной выше формулировке, фиг06, мы наблюдали значение альфа.

скажем alpha == 0, затем R’_t = R’_t-1… (мы игнорируем фактическое значение)

скажите alpha == 1, затем R’_t = R_t-1…. (мы игнорируем прогнозируемое значение)

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

Предположим, что альфа = 0,8,

R′t = α∗R’_t−1 + (1−α)∗R_t-1

R′t = (0.8 ∗ R’_t−1) + (0.2 ∗ R_t-1)

10. Базовая модель Результат

Мы ясно видим, что функция предыдущего значения работает лучше.

11. Применение регрессионных моделей машинного обучения.

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

11.1. Подготовка данных

Мы собираемся разделить данные на 70–30 частей: 70% данных для поездов и 30 для тестов. Первоначально у нас есть

# Now we computed 18 features for every data point 

# 1. cluster center latitude
# 2. cluster center longitude
# 3. day of the week 
# 4. f_t_1: number of pickups that are occurred previous t-1th 10min interval
# 5. f_t_2: number of pickups that are occurred previous t-2th 10min interval
# 6. f_t_3: number of pickups that are occurred previous t-3th 10min interval
# 7. f_t_4: number of pickups that are occurred previous t-4th 10min interval
# 8. f_t_5: number of pickups that are occurred previous t-5th 10min interval
# 9. Exponential Weighted Moving Avg
# --------------------Fourier Features --------------------------
# 10. a_1 : amplitude corresponding to 1st highest fourier transform value
# 11. f_2 : frequency corresponding to 2nd highest fourier transform value
# 12. a_2 : amplitude corresponding to 2nd highest fourier transform value
# 13. f_3 : frequency corresponding to 3rd highest fourier transform value
# 14. a_3 : amplitude corresponding to 3rd highest fourier transform value
# 15. f_4 : frequency corresponding to 4th highest fourier transform value
# 16. a_4 : amplitude corresponding to 4th highest fourier transform value
# 17. f_5 : frequency corresponding to 5th highest fourier transform value
# 18. a_5 : amplitude corresponding to 5th highest fourier transform value

Примечание: функции Фурье, упомянутые выше в разделе 8. Также обратите внимание, что мы проигнорировали компонент f_1, который является компонентом постоянного тока в преобразовании Фурье, как описано в предыдущем разделе 8.

12. Результат

Мы можем наблюдать, что простая модель, такая как «Exponential WMA», работает довольно хорошо, чем сложные модели, такие как Random Forest и XGboost. Мы хотим, чтобы стоимость MAPE была как можно меньше.

13. Прогноз Холта Винтера.

Из приведенного выше анализа становится ясно, что функция «предыдущего значения» является простой и наиболее эффективной функцией, и если данные представляют собой временные ряды, мы должны попробовать Holt Winter Forecasting, это может улучшить наши результаты.

В приведенной выше формулировке уровень может указываться как ожидаемое или прогнозируемое значение. Можно сказать тренд, предположим, что есть график XY, скорость изменения X по отношению к Y называется трендом. Это может быть восходящий или нисходящий тренд. Возрастающее или убывающее значение в серии. Если значение увеличивается, то тренд восходящий, а значение уменьшается, тогда тренд нисходящий. Сезонно, если временные ряды повторялись в течение определенного интервала времени, то можно сказать, что данные являются сезонными.

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

Приведенный ниже код используется для прогнозирования Holt Winter (тройное экспоненциально взвешенное скользящее среднее).

def initial_trend(series, slen):
    sum = 0.0
    for i in range(slen):
        sum += float(series[i+slen] - series[i]) / slen
    return sum / slen

def initial_seasonal_components(series, slen):
    seasonals = {}
    season_averages = []
    n_seasons = int(len(series)/slen)
    # compute season averages
    for j in range(n_seasons):
        season_averages.append(sum(series[slen*j:slen*j+slen])/float(slen))
    # compute initial values
    for i in range(slen):
        sum_of_vals_over_avg = 0.0
        for j in range(n_seasons):
            sum_of_vals_over_avg += series[slen*j+i]-season_averages[j]
        seasonals[i] = sum_of_vals_over_avg/n_seasons
    return seasonals


def triple_exponential_smoothing(series, slen, alpha, beta, gamma, n_preds):
    result = []
    seasonals = initial_seasonal_components(series, slen)
    for i in range(len(series)+n_preds):
        if i == 0: # initial values
            smooth = series[0]
            trend = initial_trend(series, slen)
            result.append(series[0])
            continue
        if i >= len(series): # we are forecasting
            m = i - len(series) + 1
            result.append((smooth + m*trend) + seasonals[i%slen])
        else:
            val = series[i]
            last_smooth, smooth = smooth, alpha*(val-seasonals[i%slen]) + (1-alpha)*(smooth+trend)
            trend = beta * (smooth-last_smooth) + (1-beta)*trend
            seasonals[i%slen] = gamma*(val-smooth) + (1-gamma)*seasonals[i%slen]
            result.append(smooth+trend+seasonals[i%slen])
    return result

Ключевым моментом для улучшения результата является точная настройка значений альфа, бета и гаммы. У меня есть тонко настроенные значения альфа, бета и гамма на шкале ниже.

alpha = [0.1,0.2,0.3,0.4]
beta = [0.1,0.15,0.20,0.25]
gamma = [0.1,0.3,0.4,0.5,0.65,0.75,0.85,0.95]
Best value I found for our problems is :
alpha=0.1
beta=0.1
gamma = 0.95

Используя приведенные выше значения, я добавил в список еще одну функцию. Я добавил «19-ю функцию» под названием «triple_Exp». Затем примените модель «19 функций» к Случайному лесу и XGboost.

14. Окончательные результаты

+-------------+---------------+---------+
|    Metric   | Random Forest | XgBoost |
+-------------+---------------+---------+
| Train MAPE  |    0.04707    | 0.03034 |
|  Test MAPE  |    0.04657    | 0.02932 |
+-------------+---------------+---------+

Это великолепно. Мы значительно улучшили характеристики модели. Мы ожидали, что убыток составит 0,12, и мы уменьшили убыток до 0,02, и это великолепно.

Пока что наша лучшая модель - XGboost Regressor с Test MAPE: 0,02932.

15. Заключение

В данных временных рядов лучше использовать «Moving Avg», «Weighted MA», «Exponential WMA», «Triple EWMA» в качестве функции. Также «Преобразование Фурье» может дать лучшие функции, но в нашем случае оно сработало не так, как ожидалось.

Точная настройка - это ключевой момент в повышении производительности модели.

Исходный код этого кейса можно найти здесь.

16. Ссылки

  1. Https://grisha.org/blog/2016/02/17/triple-exponential-smoothing-forecasting-part-iii/
  2. Https://scikit-learn.org/stable/modules/generated/sklearn.cluster.MiniBatchKMeans.html
  3. Https://www.youtube.com/watch?v=DUyZl-abnNM
  4. Https://xgboost.readthedocs.io/en/latest/
  5. Https://www.appliedaicourse.com/

Вы можете проверить похожие интересные блоги здесь:

Приложение для обнаружения мемов с использованием глубокого обучения.

Блог прогнозов Zomato Rate.