如何向圆形直方图添加置信区间（冯·米塞斯分布）

debugcn 发表于 Dev

科恩病毒

我有时间数据，我想在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()

结果如下图：