机器算法验证 - 平滑循环/周期性时间序列 - 吾爱随笔录

平滑循环/周期性时间序列

机器算法验证 r 密度函数内核平滑样条

2022-03-25 04:12:24

我有一天中每小时的机动车事故数据。正如您所料，它们在一天的中间很高，在高峰时段达到高峰。ggplot2 的默认 geom_density 可以很好地平滑它

与酒驾相关的车祸数据的一个子集在一天结束时（晚上和清晨）都很高，在极端情况下最高。但是 ggplot2 的默认 geom_density 仍然下降到右手边的极端。

该怎么办？目的仅仅是可视化——不需要（有吗？）进行稳健的统计分析。

伊姆古尔

x <- structure(list(hour = c(14, 1, 1, 9, 2, 11, 20, 5, 22, 13, 21, 
                        2, 22, 10, 18, 0, 2, 1, 2, 15, 20, 23, 17, 3, 3, 16, 19, 23, 
                        3, 4, 4, 22, 2, 21, 20, 1, 19, 18, 17, 23, 23, 3, 11, 4, 23, 
                        4, 7, 2, 3, 19, 2, 18, 3, 17, 1, 9, 19, 23, 9, 6, 2, 1, 23, 21, 
                        22, 22, 22, 20, 1, 21, 6, 2, 22, 23, 19, 17, 19, 3, 22, 21, 4, 
                        10, 17, 23, 3, 7, 19, 16, 2, 23, 4, 5, 1, 20, 7, 21, 19, 2, 21)
               , count = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L))
          , .Names = c("hour", "count")
          , row.names = c(8L, 9L, 10L, 29L, 33L, 48L, 51L, 55L, 69L, 72L, 97L, 108L, 113L, 
                          118L, 126L, 140L, 150L, 171L, 177L, 184L, 202L, 230L, 236L, 240L, 
                          242L, 261L, 262L, 280L, 284L, 286L, 287L, 301L, 318L, 322L, 372L, 
                          380L, 385L, 432L, 448L, 462L, 463L, 495L, 539L, 557L, 563L, 566L, 
                          570L, 577L, 599L, 605L, 609L, 615L, 617L, 624L, 663L, 673L, 679L, 
                          682L, 707L, 730L, 733L, 746L, 754L, 757L, 762L, 781L, 793L, 815L, 
                          817L, 823L, 826L, 856L, 864L, 869L, 877L, 895L, 899L, 918L, 929L, 
                          937L, 962L, 963L, 978L, 980L, 981L, 995L, 1004L, 1005L, 1007L, 
                          1008L, 1012L, 1015L, 1020L, 1027L, 1055L, 1060L, 1078L, 1079L, 
                          1084L)
          , class = "data.frame")

ggplot(x, aes(hour)) + 
  geom_bar(binwidth = 1, position = "dodge", fill = "grey") +
  geom_density() + 
  aes(y = ..count..) +
  scale_x_continuous(breaks = seq(0,24,4))

很高兴任何拥有更好统计词汇的人来编辑这个问题，尤其是标题和标签。

4个回答

我不经常使用 R 也从未使用过ggplot，但这里有一个简单的故事，或者我猜是这样。

一天中的时间显然是一个循环或周期性变量。在您的数据中，您有 0(1)23 小时环绕，因此 23 后面是 0。但是，ggplot至少从您提供的信息中不知道。就它而言，可能有 -1、-2 等或 24、25 等处的值，因此一些概率可能被平滑到超出观察数据的限制，实际上超出了可能的数据。

这也将发生在您的主要数据上，但不是很明显。

如果您想要对此类数据进行核密度估计，您需要一个足够聪明的例程来正确处理此类周期性或循环变量。“正确”意味着例程在圆形空间上进行平滑处理，认识到 0 跟在 23 之后。在某些方面，这种分布的平滑处理比通常情况更容易，因为没有边界问题（因为没有边界）。其他人应该能够就在 R 中使用的函数提供建议。

这种数据介于周期性时间序列和循环统计之间。

提供的数据有 99 个观察值。为此，直方图效果很好，尽管我可以看到您可能希望对其进行一些平滑处理。

在此处输入图像描述

（更新）这是一个品味和判断的问题，但我认为你的平滑曲线过于平滑。

这里作为样本是双权重密度估计。我使用我自己的 Stata 程序来处理循环数据，以度为单位，即席转换为 15 *（小时 + 0.5），但密度以每小时表示。相比之下，这有点不够平滑，但您可以调整您的选择。

在此处输入图像描述

要使周期性平滑（在任何平台上），只需将数据附加到自身，平滑较长的列表，并切断末端。

这是一个R插图：

y <- sqrt(table(factor(x[,"hour"], levels=0:23)))
y <- c(y,y,y)
x.mid <- 1:24; offset <- 24
plot(x.mid-1, y[x.mid+offset]^2, pch=19, xlab="Hour", ylab="Count")
y.smooth <- lowess(y, f=1/8)
lines(x.mid-1, y.smooth$y[x.mid+offset]^2, lwd=2, col="Blue")

（因为这些是我选择平滑平方根的计数；它们被转换回用于绘图的计数。）跨度lowess已从其默认值大幅缩小，f=2/3因为（a）我们现在处理的数组长了三倍，这应该导致我们将减少到，并且（b）我想要一个相当局部的平滑，以便在中间三分之一处不会出现明显的端点效应。 $f$ $2/9$

它在这些数据上做得很好。特别是，第 0 小时的异常已被平滑处理。

做 Tukey 的 4253H，在三个连接的副本上两次原始计数，然后取中间的一组平滑值，与 whuber 在计数的平方根上的低点几乎相同。
在此处输入图像描述

此外，作为更复杂的替代方案，您可能希望查看周期性样条曲线。splines您可以在 R 包a和中找到适合它们的工具mgcv。我看到的优于已经建议的方法的优点是您可以计算拟合的自由度，这对于“三份”方法并不明显。

其它你可能感兴趣的问题

上一篇偏度和峰度公式的证明/推导下一篇线性回归解释