Блог

У меня есть данные о времени, и я хочу построить график частоты в час на 24-часовых часах.

Данные преобразуются в circular , и оценки для «периодического среднего» mu и «концентрации» kappa вычисляются с mle.vonmises() помощью .

График генерируется с ggplot2 помощью , используя geom_hist() и coord_polar() . Периодическое среднее отображается на графике простым вызовом geom_vline() .

Вопрос

Я хочу нарисовать доверительный интервал в 95% вокруг среднего значения. Затем я хотел бы визуально проверить, находится ли данная временная метка (например, «22:00:00») в пределах CI или нет. Как мне это сделать с распределением фон Мизеса и ggplot2?

Приведенный ниже код показывает, как далеко я продвинулся.

Данные

 timestamps <- c("08:43:48", "09:17:52", "12:56:22", "12:27:32", "10:59:23", 
                "07:22:45", "11:13:59", "10:13:26", "10:07:01", "06:09:56", 
                "12:43:17", "07:07:35", "09:36:44", "10:45:00", "08:27:36", 
                "07:55:35", "11:32:56", "13:18:35", "11:09:51", "09:46:33", 
                "06:59:12", "10:19:36", "09:39:47", "09:39:46", "18:23:54")
  

Код

 library(lubridate)
library(circular)
library(ggplot2)

## Convert from char to hours
timestamps_hrs <- as.numeric(hms(timestamps)) / 3600

## Convert to class circular
timestamps_hrs_circ <- circular(timestamps_hrs, units = "hours", template = "clock24")

## Estimate the periodic mean and the concentration 
## from the von Mises distribution
estimates <- mle.vonmises(timestamps_hrs_circ)
periodic_mean <- estimates$mu %% 24
concentration <- estimates$kappa

## Clock plot // Circular Histogram
clock01 <- ggplot(data.frame(timestamps_hrs_circ), aes(x = timestamps_hrs_circ))  
  geom_histogram(breaks = seq(0, 24), colour = "blue", fill = "lightblue")  
  coord_polar()   
  scale_x_continuous("", limits = c(0, 24), breaks = seq(0, 24), minor_breaks = NULL)  
  theme_light()

clock01

## Add the periodic_mean
clock01   
  geom_vline(xintercept = as.numeric(periodic_mean), color = "red", linetype = 3, size = 1.25) 
  

Это дает следующий график:

Я думаю, что нашел приближение решения. Поскольку мы знаем параметры mu и kappa (соответственно. среднее периодическое и концентрация), мы знаем распределение. Это, в свою очередь, означает, что мы знаем плотности заданных временных меток, и мы можем рассчитать пороговое значение для 95% доверительного уровня.

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

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

Код

(предполагается, что код в вопросе был запущен)

 quantile <- qvonmises((1 - 0.95)/2, mu = periodic_mean, kappa = concentration)
cutoff <- dvonmises(quantile, mu = periodic_mean, kappa = concentration)

## generate a timestamp for every minute in a day
## then the transformations needed
ts_1min <- format(seq.POSIXt(as.POSIXct(Sys.Date()), 
                             as.POSIXct(Sys.Date() 1), 
                             by = "1 min"), 
                  "%H:%M:%S", tz = "GMT")
ts_1min_hrs <- as.numeric(hms(ts_1min)) / 3600
ts_1min_hrs_circ <- circular(ts_1min_hrs, units = "hours", template = "clock24")
## generate densities to compare with the cutoff
dens_1min <- dvonmises(ts_1min_hrs_circ, mu = periodic_mean, kappa = concentration)
 
## compare: vector of FALSE/TRUE
feat_1min <- dens_1min >= cutoff
df_1min_feat <- data.frame(ts = ts_1min_hrs_circ, 
                             feature = feat_1min)

## get the min and max time of the CI
CI <- df_1min_feat %>% 
  filter(feature == TRUE) %>%
  summarise(min = min(ts), max= max(ts))

CI
#   min      max
# 5.283333 14.91667
  

Используя приведенную выше информацию и используя geom_rect() ее, мы можем получить то, что хотим:

 ggplot(data.frame(timestamps_hrs_circ), aes(x = timestamps_hrs_circ))  
  coord_polar()  
  scale_x_continuous("", limits = c(0, 24), breaks = seq(0, 24), minor_breaks = NULL)  
  geom_vline(xintercept = as.numeric(CI), color = "darkgreen", linetype = 1, size = 1.5)  
  geom_rect(xmin = CI$min, xmax = CI$max, ymin = 0, ymax = 5, alpha = .5, fill = "lightgreen")  
  ggtitle(label = "Circular Histogram", subtitle = "periodic mean in red,n95%-CI in green" )  
  geom_histogram(breaks = seq(0, 24), colour = "blue", fill = "lightblue")  
  geom_vline(xintercept = as.numeric(periodic_mean), color = "red", linetype = 2, size = 1.5)  
  theme_light()
  

В результате получается следующий график:

Я надеюсь, что кто-то может извлечь из этого пользу.