r - 有没有办法让 R 中的 density() 函数使用计数与概率？

Question

有没有办法让 R 中的 density() 函数使用计数与概率？

例如，在使用直方图函数 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 的方法对我不起作用:(。

Answer 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 轴划分的精细程度(这取决于 n 到 density 的参数)。您可以使用 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"))