获取数据的概率密度

Question

我需要分析关于DSL线路的互联网会话数据。我想查看会话持续时间的分布情况。我想到一个简单的方法是首先制作所有会话持续时间的概率密度图。

我已经在R中加载了数据并使用了`density（）`函数，代码如下：

plot(density(data$duration), type = "l", col = "blue", main = "Density Plot of Duration",
     xlab = "duration(h)", ylab = "probability density")

我是新手，对R和这种分析一窍不通。通过谷歌搜索，我找到了以下内容。我得到了一个图表，但是还有一些问题。这个函数是我想要做的正确函数吗？还是有其他的函数？

在图表中，我发现Y轴刻度从0到1.5。我不明白为什么是1.5，难道不应该是从0到1吗？

此外，我想获得更平滑的曲线。由于数据集非常大，所以线条非常粗糙。在展示时将它们平滑处理会更好。我该如何做呢？

- sfactor

---

你误解了密度的概念。X的密度可以看作是与从总体中抽取一个接近X的数字的机会成比例的值。现在根据定义，密度函数的积分等于1。这并不意味着密度函数的最大值应该是1，它可以很容易地更大。事实上，对于自由度为（1,1）的F分布，密度函数的最大值（在0处）甚至可以是无穷大。 - Joris Meys

@Joris 是的，我现在意识到我没有正确理解它。我过于简单地假设，由于它是概率分布，所以它会小于1 :)。 - sfactor

@JorisMeys，当曲线下的总面积为1时，如何使概率密度函数大于1？如果概率密度函数超过1，那么很可能分布不是正常的，需要进行归一化处理。 - karthiks

@karthiks 因为例如高度为10，宽度为0.01的矩形面积为0.1，而Y值（因此PDF）将为10。对于面积，您需要考虑X轴和Y轴，而不仅仅是Y轴。 - Joris Meys

Answer 1

正如nico所说，你应该查看hist，但你也可以将它们结合起来。然后你可以使用lines调用密度。

duration <- rpois(500, 10) # For duration data I assume Poisson distributed
hist(duration,
   probability = TRUE, # In stead of frequency
   breaks = "FD",      # For more breaks than the default
   col = "darkslategray4", border = "seashell3")
lines(density(duration - 0.5),   # Add the kernel density estimate (-.5 fix for the bins)
   col = "firebrick2", lwd = 3)

你应该得到类似这样的结果：

请注意，核密度估计默认假设高斯核。但带宽通常是最重要的因素。如果你直接调用 density 函数，它会报告默认的带宽估计值：

> density(duration)

Call:
        density.default(x = duration)

Data: duration (500 obs.);      Bandwidth 'bw' = 0.7752

       x                 y            
 Min.   : 0.6745   Min.   :1.160e-05  
 1st Qu.: 7.0872   1st Qu.:1.038e-03  
 Median :13.5000   Median :1.932e-02  
 Mean   :13.5000   Mean   :3.895e-02  
 3rd Qu.:19.9128   3rd Qu.:7.521e-02  
 Max.   :26.3255   Max.   :1.164e-01  

这里是0.7752。按照nico建议，为您的数据进行检查并尝试玩弄它。您可能需要查看?bw.nrd。

- eyjo

Answer 2

您应该尝试调整带宽（bw）参数以更改曲线的平滑程度。通常情况下，R会自动提供漂亮而平滑的曲线，但也许对于您特定的数据集来说情况并非如此。

至于您正在使用的调用，是的，它是正确的，type="l"不是必需的，它是用于绘制密度对象的默认值。曲线下面积（即您的密度函数从-无穷大到+无穷大的积分）将等于1。

现在，在您的情况下使用密度曲线是最好的选择吗？也许是，也许不是...这实际上取决于您想要进行的分析类型。可能使用hist就足够了，甚至可能更具信息量，因为您可以选择特定的持续时间区间（请参见?hist获取更多信息）。

- nico

---

谢谢，我会看一下，但我仍然不明白为什么密度轴会大于1。 - sfactor

正如我所说，曲线下的面积（即sum(dx*y)）等于1。y轴的实际值取决于带宽。较小的带宽值将生成较高的y值。尝试绘制density(rnorm(1000), 0.2)和density(rnorm(1000), 2)以查看差异。 - nico

直方图相对于密度看起来是右偏的。这是因为使用正态核函数假设泊松分布变量吗？ - David LeBauer

@David：我不确定R如何计算密度估计，可能也是直方图分箱的问题，但我会把答案留给比我更有经验的人。 - nico

Answer 3

我本来想把这个作为对之前回答的评论，但是它太长了。 显然的偏斜是由于直方图中数值的分组方式。使用直方图处理离散数据通常是错误的。请参见下面...

set.seed(1001)
tmpf <- function() {
  duration <- rpois(500, 10) # For duration data I assume Poisson distributed
  hist(duration,
       probability = TRUE, # In stead of frequency
       breaks = "FD",      # For more breaks than the default
       col = "darkslategray4", border = "seashell3",
       main="",ann=FALSE,axes=FALSE,xlim=c(0,25),ylim=c(0,0.15))
  box()
  lines(density(duration),   # Add the kernel density estimate
        col = "firebrick2", lwd = 3)
  par(new=TRUE)
  plot(table(factor(duration,levels=0:25))/length(duration),
       xlim=c(0,25),ylim=c(0,0.15),col=4,ann=FALSE,axes=FALSE)
}

par(mfrow=c(3,3),mar=rep(0,4))
replicate(9,tmpf())

- Ben Bolker

---

是的，没错，箱子始终会在整数的两侧（右侧=TRUE vs. 右侧=FALSE）。我主要只是用它来预先可视化数据，没有什么大碍。但是可以通过简单的-0.5密度轻松修复... - eyjo

@eyjo：这是假设您正在使用整数断点，但您不受此限制。 - nico