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

Редкие события по определению очень редки. Но они неизбежно случаются, и когда они случаются, они имеют огромные последствия. 9/11 было последним событием. Финансовый кризис 2007/08 г. стал последним событием. Коронавирус был последним событием. Многие из продуктов, которыми вы пользуетесь каждый день, услуг, с которыми вы взаимодействуете, и компаний, в которых вы работаете, являются хвостовыми событиями. Так что да, хвостовые события редки, но когда они случаются, их влияние огромно.

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

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

Теория экстремальных ценностей

Теория экстремальных значений (EVT) - это отрасль статистики, которая имеет дело с значениями, которые значительно отклоняются от центральной тенденции распределения вероятностей. Стоит отметить, что задача с экстремальной стоимостью принципиально отличается от задач прогнозирования редких событий. Например, обнаружение мошенничества - это тип проблемы классификации редких событий, но это не проблема особого значения. Почему? Потому что такое редкое событие, как мошенничество, является редким, но предсказуемым, тогда как тип экстремальных значений, с которыми связан EVT, обычно трудно предсказать.

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

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

Распределение Гамбеля похоже на нормальное распределение в том, что оно определяется исключительно двумя параметрами - местоположением и масштабом. Но, в отличие от Нормального распределения, оно асимметрично. Это означает, что хвост распределения может отклоняться в сторону экстремальных значений. На практике это означает, что распределение Гамбеля присваивает большую вероятность большему количеству экстремальных событий (т. Е. Событий или значений в хвосте распределения), чем нормальное распределение. Вы можете увидеть это на графике ниже.

np.random.seed(123)
n = 100000
fs = {"fontsize" : 16}
plt.figure(figsize=(20, 12))
sns.kdeplot(stats.distributions.gumbel_r().rvs(n), label="Gumbel") 
sns.kdeplot(stats.distributions.norm().rvs(n), label="Normal")
plt.legend(**fs)
plt.xticks(**fs)
plt.yticks(**fs)
plt.ylabel("Density", **fs)
plt.xlabel("X", **fs);

Оценка вероятности экстремальных значений

Чтобы продемонстрировать EVT на практике, я воспользуюсь этим набором данных по скоростям ветра, который включает экстремальные значения. Данные представлены ниже. В этих данных можно увидеть сильную сезонную закономерность, а также периодические значения, которые значительно выходят за пределы нормального диапазона. Первое, что нам нужно сделать, это выделить эти значения. Для этого мы можем использовать подход максимума блока. Этот метод просто извлекает максимальное значение данных из диапазона (т. Е. Блока). В этом случае мы извлечем максимальное значение для каждого месяца из данных. Эти точки отмечены на графике красными точками.

block_maxima_month = (
    df
    .set_index("date")
    .resample("M")
    .max()
    .wind_speed
    .reset_index()
    .query("wind_speed > 5")
)
fs = {"fontsize" : 16}
plt.figure(figsize=(24, 10))
plt.scatter(df.date, df.wind_speed, alpha=.3)
plt.plot(df.date, df.wind_speed, alpha=.7, c="b")
plt.scatter(block_maxima_month.date, block_maxima_month.wind_speed, c="red", s=80)
plt.xticks(**fs)
plt.yticks(**fs)
plt.ylabel("Wind speed", **fs)
plt.xlabel("Date", **fs);

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

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

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

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

# Extract posterior distributions
ev_post = pm.sample_posterior_predictive(ev_trace, model=ev_model)["lik"].T
norm_post = pm.sample_posterior_predictive(norm_trace, model=norm_model)["lik"].T
# Plot density 
plt.figure(figsize=(20, 10))
pm.plot_dist(ev_post, label="Gumbel", textsize=20, hist_kwargs=fs)
pm.plot_dist(norm_post, label="Normal", textsize=20, )
sns.kdeplot(block_maxima_month.wind_speed, label="Actual")
plt.legend(fontsize=14)

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

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

Приведенный выше простой код выполняет итерацию по указанным значениям и использует апостериорное распределение Гамбеля для оценки вероятности увидеть большее значение. Эти вероятности превышения показаны ниже. Например, вероятность того, что скорость ветра превысит 20, составляет около 0,38.

Резюме

Прогнозирование экстремальных значений - интересная, но сложная задача. В этом посте показан общий подход, который вы можете использовать, чтобы приступить к решению этих проблем, но он может стать гораздо более сложным. Кроме того, в этом посте я не обусловил вероятность экстремального значения какой-либо ковариатой. Одним из интересных расширений этого анализа было бы увидеть, есть ли какой-либо смысл во включении дополнительных ковариат в модель (например, скорость ветра за предыдущий день или неделю или другие погодные условия, такие как влажность). Самый важный вывод - использование распределений экстремальных значений, таких как Гамбель. Когда вероятность значения или события значительно превышает условное среднее значение распределения, вам нужен подход, который учитывает правдоподобие хвостовых событий.

Спасибо за прочтение!

P.S. Весь код этого поста можно найти здесь.