Те, кто работает с временными рядами, часто сталкивается с двумя проблемами. Первая – нет полных данных. Вторая – битые данные, когда встречается много выбросов, шума и пропусков. Редко встречаются случаи, когда всё было бы идеально. И данных много, и можно легко найти нужные. Такое встретишь крайне редко или почти никогда.
Возникает вопрос - как решить эту проблему? Я нашёл решение. Давайте расскажу вам, как я решаю проблему битых данных, выбросов, пропусков. Какие я использовал методы, в чем их отличия, преимущества и какие я считаю самыми лучшими.
Начнём мы с первого метода – фильтра Хэмплея. В этой статье речь пойдёт именно о нём. Я постараюсь как можно проще рассказать о его особенностях и показать всё на наглядных примерах. Приступим.
Как работает фильтр Хэмплея
Для начала стоит понять, что такое фильтр Хэмплея. В интернете о нём вы мало что найдёте. По крайней мере, я встретил лишь скудную информацию. Хотя потратил много времени на поиски нужной информации о фильтре.
Главная цель Хэмплея – найти и заменить выбросы в заданном временном ряду. Для этого в своей основе он использует скользящее среднее с заданным окном. Для каждой итерации или окна фильтр вычисляет медиану и стандартное отклонение. Оно выражается в среднем абсолютном значении и обозначается как MAD.
Чтобы MAD стал последовательной оценкой стандартного отклонения надо умножить его на постоянный коэффициент k. Коэффициент зависит от распределения. Мы считаем, что данные подчиняются распределению Гаусса, поэтому берём коэффициент равным 1,4826.
Если значение медианы окна скользящего среднего больше чем х стандартных отклонений, то это – выброс.
Фильтр Хэмплея имеет 2 настраиваемых параметра:
· размер раздвижного окна
· количество стандартных отклонений, которые идентифицируют выброс
Для начала надо импортировать нужные библиотеки:
import matplotlib.pyplot as plt
import warnings
import pandas as pd
import numpy as npЗагрузить данные из csv файла:
df = pd.read_csv('data.csv')Распечатать данные
df.head()
Вы можете заметить, что выбросы в df уже помечены. Это делается для того, чтобы мы могли сравнить работу алгоритма с фактом.
Далее визуализируем наш df
plt.plot(df.x, df.y)
plt.scatter(df[df.outlier == 1].x, df[df.outlier == 1].y, c='r', label='outlier')
Теперь можно реализовывать фильтр Хэмпеля. Для этого используем 3 стандартных отклонения. Почему именно 3? Потому что этого с лихвой хватит для нашего временного ряда.
def hampel(y, window_size, simg=3):
n = len(y)
new_y = y.copy()
k = 1.4826
idx = []
for i in range((window_size),(n - window_size)):
r_median = np.median(y[(i - window_size):(i + window_size)]) #скользящая медиана
r_mad = np.median(np.abs(y[(i - window_size):(i + window_size)] - r_median)) #скользящий MAD
if (np.abs(y[i] - r_median) > simg * r_mad):
new_y[i] = r_median #замена выброса
idx.append(i)
return new_y, idxВызываем фильтр Хэмплея с окном скользящего среднего равного 3, чтобы определить выброс. Этого будет достаточно для нашей задачи.
new_y, outliers = hampel(df.y, 3)В переменной new_y лежит новый временный ряд без выбросов. В outliers - индексы выбросов во временном ряду.
Заливаем новый временный ряд в df вместе с признаками выбросов.
df['new_y'] = new_y
df.loc[outliers, 'outlier_hampel'] = 1Осталось визуализировать данные.
from matplotlib.pyplot import figure
figure(figsize=(15, 6), dpi=80)
plt.plot(df.x, df.y)
plt.plot(df.x, df.new_y)
plt.scatter(df[df.outlier == 1].x, df[df.outlier == 1].y, c='r', label='outlier')
plt.scatter(df[df.outlier_hampel == 1].x, df[df.outlier_hampel == 1].y, c='b', label='outlier')Выбросы, размеченные вручную, выделяются красным цветом. Синие выбросы – это определение модели.
На графике видно, что красного цвета нет. Вывод – алгоритм работает на отлично.
Для лучшего понимания возьмём другой временной ряд, чтобы снова проверить алгоритм.
Здесь заметим, что выбросов гораздо больше.
new_y, outliers = hampel(df_new.y, 3)
df_new['new_y'] = res
df_new.loc[detected_outliers, 'outlier_hampel'] = 1
from matplotlib.pyplot import figure
figure(figsize=(15, 6), dpi=80)
plt.plot(df_new.x, df_new.y)
plt.plot(df_new.x, df_new.new_y)
plt.scatter(df_new[df_new.outlier == 1].x, df_new[df_new.outlier == 1].y, c='r', label='outlier')
plt.scatter(df_new[df_new.outlier_hampel == 1].x, df_new[df_new.outlier_hampel == 1].y, c='b', label='outlier')
Увеличим окно скользящего среднего.
new_y, outliers = hampel(df_new.y, 5)
df_new['new_y'] = res
df_new.loc[detected_outliers, 'outlier_hampel'] = 1
from matplotlib.pyplot import figure
figure(figsize=(15, 6), dpi=80)
plt.plot(df_new.x, df_new.y)
plt.plot(df_new.x, df_new.new_y)
plt.scatter(df_new[df_new.outlier == 1].x, df_new[df_new.outlier == 1].y, c='r', label='outlier')
plt.scatter(df_new[df_new.outlier_hampel == 1].x, df_new[df_new.outlier_hampel == 1].y, c='b', label='outlier')
Видно, что стало гораздо лучше.
Что касается точности, то она равна 93.33333333333333 %. Я считаю, что это отличный процент.
(df_new[df_new.outlier_hampel == 1].shape[0]/df_new[df_new.outlier == 1].shape[0])*100Что в итоге?
Фильтр Хэмпеля прекрасно справляется со своей задачей. Его главным преимуществом стала простота реализации. Он может работать быстро как на малых, так и на больших объемах данных. Само собой, есть что улучшить, но в качестве простого и рабочего инструмента фильтр Хэмпеля показывает себя весьма неплохо.
Но так ли он хорош, если сравнивать его с другими алгоритмами? Например, с тестом Греббса, критерием выбора Рознера и рандомным лесом. Об этом я расскажу в следующих статьях, а в конце сравним результаты работы каждого алгоритма и вынесем окончательный вердикт, что же лучше.
Комментарии (10)
sshikov
03.12.2022 13:06+2Прогнозирование продаж Python
Вообще, это легко можно воспринять так, что вы продаете Python. Почем он у вас? И как идут продажи? Как по мне, тут не хватает слов. «Прогнозирование продаж, реализация на Python», или что-то в этом духе было бы лучше.
Фильтр Хэмпеля
фильтра Хэмплея
А вот всего два разных написания для основного предмета статьи — маловато. Давайте еще вариантов придумаем? :)
adeshere
03.12.2022 17:38А вот всего два разных написания для основного предмета статьи — маловато. Давайте еще вариантов придумаем? :)
Я бы еще предложил вариант Хампель. ;-)
Сама идея оценки выбросов в скользящем окне правильная и хорошая. На эту тему существуют десятки подходов и методов. Их основная идея состоит в том, что критерий идентификации выбросов адаптируется к текущим свойствам сигнала. Чаще всего текущее значение ряда сравнивается со средним или с медианным значением в окне, а величина выброса измеряется в стандартных отклонениях, рассчитанных в том же окне либо по всему ряду в целом. Мы обычно применяем именно второй вариант, так как у нас в геофизических временных рядах выбросы часто идут пачками и в окно иногда попадает сразу несколько выбросов, что приводит к росту локальной (скользящей) дисперсии. А слишком сильная изменчивость уровня отсечения выбросов - это не всегда хорошо.
Но иногда более подходящим может оказаться и критерий скользящей дисперсии. Для многих сигналов, в том числе и в эконометрике, разброс пропорционален текущему среднему. Тогда дисперсию лучше оценивать в скользящем окне, чтобы критерий сохранял однородность (в относительных единицах) на всем протяжении ряда вне зависимости от его текущего уровня. Но тут есть один важный нюанс. Если в сигнале встречаются участки быстрого роста или быстрого спада сигнала, то стандартное отклонение в таких окнах может оказаться намного больше, чем фактический разброс данных относительно линии тренда. Поэтому для аккуратного удаления выбросов лучше сначала удалить из сигнала функциональный либо непараметрический тренд (можно взять ту же скользящую медиану), и только после этого применять критерий выделения выбросов. Идентифицированные таким способом выбросы потом выбраковываются, естественно, из исходного ряда с трендом. Подробнее такой подход к фильтрации выбросов описан вот в этой статье.
Ну и последнее замечание: все методы этого типа прекрасно работают при наличии в сигнале пропусков наблюдений. Поэтому при кодировании таких процедур нет никакого смысла заранее ограничивать себя требованием полноты данных. Ведь во всех реальных рядах неизбежно присутствуют дефекты разного рода, и очень часто методологически правильнее просто исключать такие значения из анализа вместо того, чтобы заменять их какими-то искусственными (вымышленными) значениями.
Точнее, такие данные надо не исключать, а...
иногда исключать из некоторых расчетов. А в других случаях они, наоборот, могут оказаться самыми ценными наблюдениями - у нас в геофизике такое довольно часто бывает. Поэтому мы не выбраковываем в исходном сигнале все аномальные точки автоматически, если только нет твердой уверенности, что это технический брак. Вместо этого в базе данных хранится исходный сигнал, а особые точки выбраковываются (заменяются пропусками) только при проведении некоторых расчетов. А в других случаях они обрабатываются как-то иначе. То есть, определения понятий "выброс", пропуск" и "аномалия" не фиксированы, а варьируются в зависимости от задачи, которая решается в данный момент.
Поэтому при программировании фильтров, подобных описанному в статье, лучше сразу же предусмотреть возможность наличия в данных пропущенных наблюдений, чтобы их не требовалось чем-нибудь заполнять. Например, можно закодировать пропуск каким-то особым числом, которое никогда не встречается в данных. Затем во всех функциях для оценки стандартного отклонения, медианы и прочего добавляется один оператор, который проверяет каждое значение ряда. Если оно равно пропуску (функция isNan() возвращает "увы, да"), то такое значение игнорируется. Мы обычно при этом еще задаем максимальный разрешенный процент пропущенных наблюдений в окне - если фактическое количество данных в пределах некоторого окна слишком маленькое, то такое окно бракуется целиком. Добавить такую обработку в программу совсем не сложно, но она кардинально повышает гибкость методики при процессинге данных с дефектами.
MGS62
03.12.2022 13:45+1В функции def hampel забыли домножить на k = 1.4826
KainoRhine Автор
03.12.2022 18:07Да увидел
Там такая строка будет: r_mad = k * np.median(np.abs(y[(i - window_size):(i + window_size)] - r_median))
zartdinov
Возможно эти выбросы являются неплохими индикаторами)
KainoRhine Автор
Да, совершенно верно. Однако, в тех данных с которыми работаем мы, продажи на очень большие суммы (либо на очень маленькие) являются выбросами.
Например, товар X в неделю в среднем отгружается на 150 т.р. а сегодня он ушёл со склада на сумму 1 500 000 рублей. Это явно выброс. Так как такие отгрузки, в силу специфичности бизнеса, очень редки.
Более подробно расскажу в статье, где мы на рабочих (реальных) данных будем строить модель прогнозирования продаж
materiatura
Возможно, десятикратное превышение среднего и выброс. Но как определить границу, как различать выбросы при меньших отклонениях, те которые "не явные"?
adeshere
Но как определить границу, как различать выбросы при меньших отклонениях, те которые "не явные"?
А вот это уже вопрос не к алгоритму, а к эксперту. Именно по этой причине человека пока что нельзя заменить системами на основе ML, которые прекрасно находят выбросы, если критерий задан, но совершенно не приспособлены для ответов на вопросы типа "А для чего мы хотим найти и удалить эти выбросы? Что является конечной целью такой обработки?", без ответов на которые невозможно подобрать оптимальный критерий и алгоритм дефектовки.
В программе может быть хоть сто процедур для анализа и коррекции выбросов, но эта тьма инструментов никак не поможет, если за штурвалом сидит "оператор", а не эксперт. Пока у него нет содержательного понимания - что должно считаться выбросом в данном конкретном случае - нет никакого смысла формализовать алгоритм.
Другой вопрос, как к этому пониманию прийти. Ну так для этого мы и сидим у компьютера ;-) В геофизике достаточно часто срабатывает подход, основанный на численном моделировании: мы обрабатываем сигнал много раз, варьируя критерии определения выбросов и т.д., а затем вычисляем какие-то целевые параметры и оцениваем, насколько хорош результат. Но использовать этот подход надо весьма аккуратно, так как при этом очень легко можно впасть в
иллюзию значимости
Что-то не нашел сейчас ссылку на описание этого эффекта, который заключается в том, что если мы проводим многократные повторные расчеты с разными наборами случайных данных, то просто по статистике иногда будем получать значимый результат. Если в этот момент крикнуть Эврика! и сразу же бежать
делать публикацию в Natureза повышением оклада к начальнику, забыв про предыдущие неудачные опыты, то можно, мягко говоря, очень жестко сесть в лужу. А ведь у этого эффекта даже название вроде бы есть, и где-то я его даже видел... (чертов склероз!).Так что вопрос совершенно правильный... но простого ответа на него нет
materiatura
Вы хорошо и развернуто написали. Автору стоило бы ответить на эти, да и другие вопросы в статье. Хорошая статья про конкретный фильтр. Но зачем он в прогнозировании продаж? Почему продажи - это временной ряд? Не то, чтобы я против, ну, временной так временной, но почему и зачем? Почему он подчиняется закону Гаусса? Почему 10-кратная продажа выброс? Опять, я не против, хотите считать это "выбросом", - считайте. Но объясните, плиз, что значит в данном случае слово "выброс", лучше бы его не было или все-таки это неплохо, продать в 10 раз больше чем обычно? Я понимаю, что это вопросы не к Вам. Наверное, я занудствую. У меня в первом комментарии после "те" запятая пропущена.