Анализ распределения интервалов между покупками на R

Прежде чем перейти к статье, хочу вам представить, экономическую онлайн игру Brave Knights, в которой вы можете играть и зарабатывать. Регистируйтесь, играйте и зарабатывайте!

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

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

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

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

В статье я представляю:

Код на R для анализа любых временных интервалов.
Подбор экспоненциального и степенного распределения под данные с помощью метода максимального правдоподобия (MLE). Для экспоненциального я использую fitdistr() из пакета MASS, а для степенного fit_power_law() из пакета igraph.
Проверку данных на соответствие подобранному распределению с помощью теста Колмогорова-Смирнова. Я использую функцию ks.test() из пакета stats.

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

Выгрузка данных

Для начала стоит выгрузить данные из банковского приложения. Как это сделать в Тинькоф или Рокетбанк, можно подсмотреть в статье на Хабре про статистику расходов по МСС.

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

library(data.table) # для работы с дата-фреймами

tr <- fread('operations.csv')
tr_buy <- tr[`Сумма операции` < 0]

# 1. Переводим строку в числовой формат с помощью strptime()
tr_buy[, time := strptime(`Дата операции`, "%d.%m.%Y %H:%M:%S")]

# 2. Добавим в датафрейм разность между временем двух последовательных покупок (в секундах)
tr_buy <- cbind(tr_buy, as.numeric(diff(tr_buy$time)))
colnames(tr_buy)[17] <- "intervals" #Переименуем колонку для красоты

# 3. Переведем интервалы в часы и сделаем интервалы положительными
tr_buy$intervals <- tr_buy$intervals / (60*60)*(-1)
nrow(tr_buy[intervals> 24])

# 4. Оставим интервалы меньше 24 часов
tr_buy <- tr_buy[intervals < 24, ,]

Смотрим на данные глобально

После выгрузки посчитаем основные описательные статистики среди интервалов до 24 часов. Таких оказалось 1098, а 48 интервалов в рассматриваемый промежуток не попали.

Медиана (ч): 1.07
Среднее (ч): 4.03

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

Раскройте для просмотра кода на R

library(ggplot2)     # для построения графиков
library(gridExtra)   # для соединения графиков

Медиана и среднее
med <- median(tr_buy$intervals)
men <- mean(tr_buy$intervals)
1. Строим график функции плотности вероятности
g1 <- ggplot(data = tr_buy, aes(x = intervals)) + 
  geom_density(alpha = 0.5) + 
  coord_cartesian(xlim = c(0, 24)) +
  labs(x = "Time intervals between two consecutive transactions (hours)", 
       y = NULL) +
  theme_bw()  
2. Строим гистограмму распределения
g2 <- ggplot(data = tr_buy, aes(x = intervals)) + 
  geom_histogram(alpha = 0.5, binwidth = 1, colour = "black") + 
  coord_cartesian(xlim = c(0, 24)) +
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  annotate("label", x = 2.7, y = 240, 
           size = 3, 
           label= "n = 1098",
           fill="grey", alpha = 0.5) + 
  annotate("label", x = 23, y = 200, 
           size = 3, 
           label= "n = 48", 
           fill="grey", alpha = 0.5) + 
  annotate("text", x = 23, y = 120, 
           size = 3, 
           label = "(out of range)") + 
  theme_bw() 
3. Строим боксплот
g3 <- ggplot(data = tr_buy, aes(x = intervals)) + 
  geom_boxplot(alpha = 0.5)+ 
  coord_cartesian(xlim = c(0, 24))+
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  annotate("label", x = med + 1.725, y = 0.15, size = 3, label= "median = 1.07", 
           fill="grey", alpha = 0.5) + 
    annotate("label", x = med + 1.725, y = -0.15, size = 3, label= "mean = 4.03", 
           fill="grey", alpha = 0.5) + 
  theme_bw() 
4. Соединяем графики вместе
grid.arrange(g1, g2, g3, nrow = 3, 
             top = paste(
               "The distribution of intervals between two consecutive transactions of Vladimir Silkin (me)\n", "Global view"))

После анализа полученного распределения я принимаю решение разделить данные на 2 группы -- примерно от 0 до 8 часов и от 8 до 24.

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

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

Берем фокус на двух частях распределения

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

Раскройте для просмотра кода на R

# Добавим в данные качественную переменную, разделяющую наши данные на 2 распределения
tr_buy$distr <- ifelse(tr_buy$intervals < 8, 'exp', 'norm')

# Медиана и среднее для экспоненциального распределения
med_exp <- median(tr_buy[intervals < 8, intervals,])
mean_exp <- (tr_buy[intervals < 8, intervals,])

# Медиана, среднее и стандартное отклонение для нормального распределения
med_norm <- median(tr_buy[intervals > 8, intervals,])
mean_norm <- mean(tr_buy[intervals > 8, intervals,])
sd_norm <- sd(tr_buy[intervals > 8, intervals,])

# Гистограмма с цветовым разделением данных.
g1 <- ggplot(data = tr_buy, aes(x = intervals, fill = distr)) + 
  geom_histogram(alpha = 0.5, binwidth = 1, show.legend = FALSE, colour = "black") + 
  coord_cartesian(xlim = c(0, 24)) +
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  annotate("label", x = 4, y = 150, 
           size = 4, 
           label= "n = 894",
           fill="red", alpha = 0.25) + 
  annotate("label", x = 16, y = 100, 
           size = 4, 
           label= "n = 204", 
           fill="#00AFBB", alpha = 0.25) + 
  theme_bw() 

# Двойной боксплот так же с разделением.
g2 <- ggplot(data = tr_buy, aes(x = intervals, fill = distr)) + 
  geom_boxplot(alpha = 0.5, show.legend = FALSE)+ 
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  
  annotate("label", x = 15.7, y = -0.1, 
           size = 4, 
           label= "median = 15.69",
           fill="#00AFBB", alpha = 0.25) + 
  annotate("label", x = 15.7, y = -0.2, 
           size = 4, 
           label= "mean = 15.77", 
           fill="#00AFBB", alpha = 0.25) + 
  annotate("label", x = 15.7, y = -0.3, 
           size = 4, 
           label= "sd = 3.88", 
           fill="#00AFBB", alpha = 0.25) + 
    
  annotate("label", x = 2, y = 0.25, 
           size = 4, 
           label= "median = 0.53",
           fill="red", alpha = 0.2) + 
  annotate("label", x = 2, y = 0.15, 
           size = 4, 
           label= "mean = 1.35", 
           fill="red", alpha = 0.2) + 
  

  
  coord_cartesian(xlim = c(0, 24)) +
  theme_bw() 

# 4. Соединяем график вместе
grid.arrange(g1, g2, nrow = 2, top = "The distribution of intervals between two consecutive transactions of Vladimir Silkin (me)\n Focus on two parts of distribution")

1. Проверка на нормальность

Сперва разберемся с околонормальным распределением. Посчитаем вероятность получить такие как у нас и более выраженные отклонения от нормальности с помощью тета Шапиро-Уилка.

shapiro.test(tr_buy[intervals > 8, intervals])

Получаем $p.value \approx 0.01$ То есть, если наша выборка, действительно, взята из нормального распределения, то вероятность получить такие и сколько угодно более выраженные отклонения от нормальности случайно равняется 1%. Я не стану отклонять гипотезу о том, что наши данные взяты из генеральной совокупности с нормальным распределением признака.

Раскройте для просмотра кода на R

g1 <- ggplot(data = tr_buy[intervals > 8], aes(x = intervals)) + 
  geom_histogram(aes(y = stat(count) / sum(count)), fill = "#00AFBB",
                 alpha = 0.4, 
                 colour = "black",
                 binwidth = 1, 
                 show.legend = FALSE) + 
  geom_density(size = 1) +
  coord_cartesian(xlim = c(8, 24)) +
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  stat_function(fun = dnorm,
                args = list(16,4), 
                size = 1, 
                colour = 'blue',
                lty = "dashed") +
  
  theme_bw() 

grid.arrange(g1, nrow = 1, top = "The distribution of intervals between two consecutive transactions of Vladimir Silkin (me) \n Long intervals (longer than 8 hours)")

На графике сплошная линия показывает функцию плотности наших данных, а пунктирная - функцию плотности идеального нормального распределения.

2. Подбор экспоненциального и степенного распределений.

Теперь попробуем фитануть плотность вероятности с помощью двух распределений - экспоненциального и степенного. Они перед вами.

$\begin{gather} f_{exp}(x) = \alpha \, e^{-\alpha \, x} \\ f_{pow}(x) = x^{-\alpha} \end{gather}$

В работе Data-Scientist`а нет ничего сложного. Главное -- подгрузить нужные библиотеки. (Маэстро)

library(MASS) # Для подбора экспоненциального распределения с помощью fitdistr()
library(igraph) # Для подбора степенного распределения с помощью fit_power_law()
library(stats) # Для проведения теста Колмогорова-Смирнова

Фитуем данные экспоненциальным распределением:

fit_e <- fitdistr(tr_buy[intervals < 8, intervals], "exponential") 
fit_e$estimate[1][[1]] # alpha
fit_e$sd[1][[1]] # alpha_error
ks.test(tr_buy[intervals < 8, intervals], "pexp", fit_e$estimate)$p.value #p.value

Получим $\alpha \approx 0.74$ и $p.value \approx 10^{-14}$

Первое найдено с помощью метода максимального правдоподобия, а второе - с помощью теста Колмогорова-Смирнова.

Физический смысл результатов теста такой же, как и у результатов теста Шапиро-Уилка: если в генеральной совокупности данные, действительно, распределены экспоненциально, то вероятность получить такие и сколь угодно выраженные отклонения от экспоненциального распределения случайно, практически равна нулю. Про рассчитываемую в тесте статистику можно почитать в англоязычной статье википедии.

Фитуем данные степенным распределением:

fit_p <- fit_power_law(tr_buy[intervals < 8, intervals], implementation = 'plfit')
fit_p$alpha # alpha
fit_p$KS.p  # p-уровень значимости из теста Колмогорова-Смирнова

Получим $\alpha \approx 2.17$ и $p.value \approx 10^{-5}$

В случае со степенным распределением p.value повышается, но все равно не достигает ни одного общепринятого порога значимости.

Составим итоговую таблицу с результатами:

Equation	Parameter	Significance
$f_{exp}(x) = \alpha \, e^{-\alpha \, x}$	$\alpha \approx 0.74$	$p.value \approx 10^{-14}$
$f_{pow}(x) = x^{-\alpha}$	$\alpha \approx 2.17$	$p.value \approx 10^{-5}$

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

Раскройте для просмотра кода на R

plaw <- function(x, a) {x^(-a)}

gf <- ggplot(data = tr_buy[intervals < 8], aes(x = intervals)) + 
  geom_density(size = 1) + 
  geom_histogram(alpha = 0.3, 
                 aes(y = stat(count) / sum(count)), 
                 binwidth = 0.5,
                 fill = "red",
                 colour = "black")+
  
  stat_function(fun = dexp, args = list(0.738317), size = 1, colour = 'blue') +
  stat_function(fun = plaw, args = list(2.171511), size = 1, colour = 'red') +

  
  coord_cartesian(xlim = c(0, 8),
                  ylim = c(0, 0.5))+
  labs(x = "Time intervals between two consecutive transactions (hours)", y = NULL) +
  theme_bw() 

grid.arrange(gf, nrow = 1, top = "The distribution of time between two consecutive transactions of Vladimir Silkin (me) \n Short intervals (less than 8 hours)")

Выводы

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

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

Есть ли какой-то физический смысл у медианы распределения интервалов? Сейчас мне кажется, что она характеризует способность человека принимать спонтанные решения о покупках, поэтому я ожидаю, что она будет увеличиваться с возрастом.
Можно ли что-то сказать о теоретическом распределении интервалов между транзакциями и можно ли на основе полученных результатов утверждать о принадлежности к нему наших данных?

Источник: https://habr.com/ru/post/575116/

Вернуться к списку

Интересные статьи

Линейный криптоанализ на примере блочного алгоритма шифрования NUSH

В современном мире остро стоит вопрос о конфиденциальности данных при обмене ими и их хранении, которая достигает за счет все возможных способов шифрования. Однако при по...

Квантовая криптография: простейшие протоколы и чуть-чуть криптоанализа

Самая известная криптографическая проблема - передача секретных сообщений. Для этой задачи чаще всего используют криптосистемы с закрытым ключом: Алиса (отправит...

Оптимизация распределения серверов по стойкам

В одном из чатов мне задали вопрос: — А есть что-то почитать, как правильно упаковывать сервера в стойки? Я понял, что такого текста не знаю, поэтому написал свой. Во-первых, этот текст...

Ищем и анализируем ошибки в коде Orchard CMS

Эта статья – результат повторной проверки проекта Orchard с помощью статического анализатора PVS-Studio. Orchard это система управления контентом с открытым исходным кодом, которая является час...

Southbridge в Челябинске и Битрикс в Kubernetes

В Челябинске проходят митапы системных администраторов Sysadminka, и на последнем из них я делал доклад о нашем решении для работы приложений на 1С-Битрикс в Kubernetes. Битрикс, Kubernetes, Сep...