有没有办法让 R 中的密度（）函数使用计数与概率？

Posted 2023-02-14

【中文标题】有没有办法让 R 中的密度（）函数使用计数与概率？

【英文标题】：Is there a way to make the density() function in R use counts vs. probability?

【发布时间】：2015-12-08 14:22:03

【问题描述】：

例如，在使用直方图函数hist 检查密度分布时，我有两个选择：

hist(x,freq=F)  #"graphic is a representation of frequencies, the counts component of the result"
hist(x,freq=T)  #"probability densities, component density, are plotted (so that the histogram has a total area of one)"

我想知道是否有办法使用density 函数做类似的事情？

在我的具体示例中，我计算了不同直径的树木。 （我会注意到，我将数据保留为连续的大小比例，而不是将它们集中到离散的大小类别中）。当我将density 函数与此数据（即plot(density(dat$D,na.rm=T,from=0))）一起使用时，它为我提供了每种尺寸概率的密度估计（当然是平滑的）。我更有兴趣将这些数据报告为茎/面积与概率，所以我更喜欢密度估计来使用计数。

想法？？

---

更新：

以下是一些真实的示例数据：

 dat <- c(6.6, 7.1, 8.4, 27.4, 11.9, 18.8, 8.9, 25.4, 8.9, 8.6, 11.4, 19.3, 7.6, 42.2, 20.8, 25.1, 38.1, 42.2, 5.2, 34.3, 42.7, 34, 37.3, 45.5, 39.4, 25.1, 30.7, 23.1, 43.4, 19.6, 30.5, 23.9, 10.7, 18.3, 30, 35.8, 8.1, 11.9, 28.4, 30.5, 34.3, 10.4, 45, 38.9, 8.9, 11.7, 9.7, 7.4, 3.8, 20.6, 48.8, 6.6, 40.4, 13, 16, 8.6, 16, 13, 12.2, 11.4, 10.2, 22.6, 17.3, 12.4, 9.7, 17.3, 10.9, 27.2, 9.1, 13, 10.9, 15, 10.4, 27.2, 21.6, 18.8, 12.7, 15.5, 17, 16.3, 18, 26.9, 10.2, 21.3, 19, 11.7, 10.7, 18, 9.9, 16.5, 19.6, 22.1, 9.9, 18.3, 17, 6.9, 7.6, 12.7, 13.2, 9.7, 13.5, 18.3, 19.3, 30, 20.1, 18.5, 12.2, 16, 17, 14.2, 5.6, 12.2, 7.6, 17, 14, 16.5, 13.7, 11.9, 14.2, 15, 13.7, 13.2, 9.1, 6.9, 9.9, 11.4, 12.7, 10.2, 12.4, 15, 20.1, 6.9, 8.1, 11.4, 10.7, 10.9, 18.3, 9.1, 6.3, 17.3, 20.1, 9.4, 7.1, 16, 15, 10.9, 14.7, 18.8, 14.5, 10.7, 14, 10.4, 14.5, 15.7, 10.9, 14.7, 19.3, 12.4, 7.1, 14, 15.5, 36.8, 23.1, 7.9, 9.9, 8.1, 14.7, 13.7, 18, 10.7, 11.9, 12.7, 12.4, 17.8, 7.9, 12.2, 10.4, 13, 14.7, 12.7, 8.1, 14.2, 10.2, 11.9, 5.6, 8.4, 6.1, 7.6, 7.9, 19.8, 7.4, 12.7, 10.2, 12.4, 10.4, 12.4, 26.9, 12.7, 16.8, 22.9, 15.7, 10.4, 13.7, 8.1, 13.7, 14.2, 21.6, 20.8, 12.4, 10.9, 10.2, 29.5, 19.3, 8.9, 6.1, 11.2, 7.1, 28.7, 15.7, 10.4, 8.6, 10.4, 9.1, 14.5, 25.7, 11.4, 15.5, 8.1, 13.2, 16.8, 5.8, 20.8, 10.2, 9.1, 5.6, 14.5, 14.5, 17.5, 29.2, 13, 14, 12.4, 9.9, 21.1, 18.8, 14, 15.5, 9.7, 24.1, 20.1, 20.3, 12.4, 15.2, 15.7, 8.6, 8.6, 10.4, 12.4, 16.8, 4.1, 8.1, 6.6, 11.7, 7.9, 17.5, 9.1, 4.6, 7.1, 7.6, 9.4, 20.8, 11.4, 15.5, 7.1, 18.5, 7.9, 16.5, 6.3, 6.1, 16.5, 15.5, 17.3, 20.3, 12.7, 20.3, 13.7, 8.4, 16.8, 14, 18, 10.9, 19.8, 10.7, 27.2, 11.4, 7.9, 11.2, 14.5, 14.2, 11.2, 13.5, 18.5, 4.3, 7.9, 6.1, 9.9, 14.7, 8.4, 14, 12.4, 15, 14.2, 11.4, 7.6, 12.7, 5.8, 16, 7.9, 3.3, 5.8, 4.8, 4.8, 7.4, 9.1, 8.4, 3.8, 9.1, 9.4, 8.4, 9.9, 7.9, 13.2, 20.8, 18.3, 16.8, 13.5, 12.4, 8.1, 6.3, 7.6, 18.5, 14, 10.2, 9.4, 11.9, 11.4, 13, 14.5, 17, 7.9, 10.2, 7.4, 5.3, 6.9, 17.8, 5.6, 10.9, 9.9, 9.9, 16.5, 8.9, 24.1, 22.9, 13.5, 10.7, 23.4, 10.9, 28.2, 5.6, 19.6, 15.2, 6.3, 23.1, 19.3, 26.7, 30.5, 13.7, 7.9, 20.8, 19.8, 21.6, 21.6, 9.9, 30.5, 16.3, 11.9, 5.1, 15.2, 13.2, 7.1, 5.8, 9.9, 19.3, 15.5, 25.7, 14, 29.7, 11.9, 12.7, 25.9, 16.3, 25.9, 6.1, 26.7, 7.9, 9.7, 22.1, 20.1, 24.4, 17.3, 13.2, 16.5, 16.8, 21.8, 15.2, 9.9, 19.6, 23.6, 23.4, 17.8, 15.5, 11.4, 20.8, 22.1, 26.4, 12.4, 14.2, 6.9, 22.1, 22.6, 34.5, 15, 13.2, 19.6, 18.3, 15.5, 13.5, 14, 19.8, 21.1, 16.3, 19.8, 13.7, 12.2, 11.7, 31.7, 12.7, 13.2, 7.6, 12.2, 13.2, 31.7, 9.9, 10.2, 9.1, 9.1, 21.6, 8.6, 12.7, 13.5, 9.7, 8.9, 11.7, 8.4, 19.6, 7.6, 13.2, 18.3, 11.2, 22.4, 10.9, 14.7, 12.7, 16.8, 18.8, 15, 8.1, 20.8, 22.1, 7.6, 16.3, 10.9, 8.9, 11.7, 24.4, 29, 29.2, 27.4, 25.1, 6.6, 11.7, 16.5)

这里正在尝试@eipi10 suggests的方法：

#Produce graph showing counts of values using table():
  plot(x=names(table(dat)), y = table(dat),type='l')
#Produce graph showing counts of values using density + @eipi10's method
  dens <- density(x = dat, na.rm = T, bw = 0.1, n = length(dat))
  dens$y <- length(dat)/sum(dens$y) * dens$y  #"fix" to counts
  plot(dens)

此代码创建以下 2 个图表 [标题为 post-hoc]：

如您所见，这两种方法在 y 轴上产生了不同的值。换句话说，@eipi10 的方法对我不起作用 :(.

【问题讨论】：

如果将它们分成离散的类，那么您所要做的就是整合每个类的曲线。

这是一个非常接近的匹配：plot(x=as.numeric(names(table(dat))), y = table(dat),type='l'); lines(dens$x,dens$y*sum(dens$y)/diff(dens$x)[1],col=2)

@BenBolker 感谢您提供替代方法。但是，我无法让它与其他数据集保持一致。无论如何，最终，我不只是希望这两个图表看起来相似。我希望 实际上 能够将 density() 输出的概率值转换为 actual 计数。这有可能吗？

【参考方案1】：

您可以通过将密度值标准化为样本中值的数量来转换为计数。例如：

# Fake data
k=1000
set.seed(104)
val = rnorm(k)
dens = density(val, n=512)

# Convert to counts
dens$y = k/sum(dens$y) * dens$y

plot(dens)

但请记住，您最终得到的计数取决于您划分 x 轴的精细程度（这取决于 density 的 n 参数）。您可以使用mean(diff(dens$x)) 确定 delta-x（间隔并没有真正变化，但由于舍入误差，它们并不完全相同）。

更新：根据您的评论，下面的代码应该可以解释发生了什么。但首先，请注意，在对实际数据进行分箱时获得的计数（通常）与从核密度估计得出的计数不匹配，除非实际数据的分箱间隔与用于核密度估计的间隔相同。 （由于核密度估计的平滑，计数在任何情况下都不太可能完全匹配，但分箱间隔需要相同才能获得密切对应。）

library(ggplot2)
library(reshape2)
library(dplyr)

# Fake data
k=1000
set.seed(104)
dat = data.frame(diameter = rnorm(k,100,10))

创建 3 个核密度估计值：前两个分别使用 20 和 100 个点。第三个使用 100 个点，但默认带宽为 1/10。

# Convert density to counts
ctc = function(data, nPoints, numValues, adj=1) 
  dens = density(data$diameter, n=nPoints, adjust=adj)
  dens$y = numValues/sum(dens$y) * dens$y
  return(dens)


dens20 = ctc(dat, 20, k)
dens100 = ctc(dat, 100, k)
dens100adj = ctc(dat, 100, k, 0.1)

使用实际计数和根据内核密度估计估计的计数创建数据框。我们将使用cut 函数来确保实际计数使用与核密度估计值相同的间隔。

dd = function(data, dens) 
  data = data.frame(table(cut(data$diameter, 
                              breaks=c(dens$x - 0.5*mean(diff(dens$x)),Inf))),
                    DensityCounts=round(dens$y,1))  # Rounding is just for easier comparison by eye if you display the data frame
  names(data)[1:2] = c("DiameterRange","ActualCounts")
  return(data)


dat20 = dd(dat, dens20)
dat100 = dd(dat, dens100)
dat100adj = dd(dat, dens100adj)

现在我们创建将每个核密度估计值与实际计数进行比较的图。请注意实际计数何时与根据密度估计创建的计数相匹配，以及带宽如何影响我们使用的间隔。

pf = function(data, title) 
  ggplot(data %>% melt(id.var="DiameterRange"), 
         aes(DiameterRange, value, colour=variable, group=variable)) +
    geom_line() +
    theme(axis.text.x=element_text(angle=-90, vjust=0.5, hjust=0)) +
    ggtitle(title)


gridExtra::grid.arrange(pf(dat20, "n=20"), 
                        pf(dat100, "n=100"), 
                        pf(dat100adj, "n=100; 1/10th default bandwidth"))

【讨论】：

这似乎不起作用。当我计数为 7 时，该 x 值的 dens$y 为 0.2865707498，但当我使用您的方程时，结果为 2.544275344。我希望结果是 7。我只是错误地接近这个吗？？

如果您发布数据样本可能会有所帮助。当您说“计数为 7”时，在什么间隔内（在您的情况下为直径范围）？间隔与您的密度估计中的间隔相同吗？如果不是，则计数可能不同。实际上，由于核密度估计的平滑，计数在任何情况下都会有所不同。

请记住，密度估计并不是 X 值（在这种情况下 X 是直径）恰好等于某个值（在这种情况下为 0.2865）的概率。这是它落在 Xi 和 Xi+1 之间的概率，其中该范围是（在核密度估计的情况下）mean(diff(dens$x))。当我们转换为计数时，我们会得到给定计数落在 Xi 和 Xi+1 范围内的概率（给定样本中的数据点数量及其特定值）。

（带宽我们没讲过，但这也会影响结果。R用默认方法确定，但你可以改变它。）

嘿@eipi10，我终于决定重新审视这个。我已经对我的问题进行了更新，包括提供示例数据和演示我如何无法使用您的方法来处理我的数据。非常感谢您提供的任何其他帮助！

【参考方案2】：

除非你特别需要density函数，否则你可以使用函数table

Counts<-table(factor(dat$D,levels=0:n)) # n=number of size levels

plot(Counts,type="l")

这将为您提供每个直径的计数，但不会平滑。

【讨论】：

我实际上是在使用 table(dat) 来通知我的图表。但是，我确实希望线条平滑