具有总和为1的密度曲线的ggplot2直方图[关闭]

Question

为非标准化数据绘制密度曲线总和为 1 的直方图非常困难。可笑。对此已经有很多问题，但他们的解决方案都不适用于我的数据。需要有一个简单有效的解决方案。我找不到简单有效的解决方案。

一些例子：

解决方案仅适用于标准化的正常数据 ggplot2: Overlay histogram with density curve

有离散数据，没有密度曲线 ggplot2 density histogram with width=.5, vline and centered bar positions

没有答案 Overlay density and histogram plot with ggplot2 using custom bins

我的数据中的密度总和不等于 1 Creating a density histogram in ggplot2?

我的数据总和不等于 1 ggplot2 density histogram with custom bin edges

这里有例子的长解释，但我的数据密度不是 1 "Density" curve overlay on histogram where vertical axis is frequency (aka count) or relative frequency?

--

一些示例代码：

#Example code
set.seed(1)
t = data.frame(r = runif(100))

#first we try the obvious simple solution that should work
ggplot(t, aes(r)) + 
  geom_histogram() + 
  geom_density()

所以，显然密度总和不等于 1。

#maybe geom_histogram needs a ..density.. ?
ggplot(t, aes(r)) + 
  geom_histogram(aes(y = ..density..)) + 
  geom_density()

它确实改变了一些东西，但不正确。

#maybe geom_density needs a ..density.. too ?
ggplot(t, aes(r)) + 
  geom_histogram(aes(y = ..density..)) + 
  geom_density(aes(y = ..density..))

那里没有变化。

#maybe binwidth = 1?
ggplot(t, aes(r)) + 
  geom_histogram(aes(y = ..density..), binwidth=1) + 
  geom_density(aes(y = ..density..))

仍然是错误的密度曲线，但现在直方图也是错误的。

可以肯定的是，我确实花了 4 个小时尝试了 ..count.. 和 ..sum.. 和 ..density.. 的各种组合，但是因为我找不到任何关于这些是如何假设的文档要工作，这是半盲试错。

所以我放弃了，避免使用ggplot2来汇总数据。

所以首先我们需要得到正确的比例data.frame，这不是那么简单：

get_prop_table = function(x, breaks_=20){
  library(magrittr)
  library(plyr)
  x_prop_table = cut(x, 20) %>% table(.) %>% prop.table %>% data.frame
  colnames(x_prop_table) = c("interval", "density")
  intervals = x_prop_table$interval %>% as.character
  fetch_numbers = str_extract_all(intervals, "\\d\\.\\d*")
  x_prop_table$means = laply(fetch_numbers, function(x) {
    x %>% as.numeric %>% mean
  })
  return(x_prop_table)
}

t_df = get_prop_table(t$r)

这给出了我们想要的那种汇总数据：

> head(t_df)
          interval density    means
1 (0.00859,0.0585]    0.06 0.033545
2   (0.0585,0.107]    0.09 0.082750
3    (0.107,0.156]    0.07 0.131500
4    (0.156,0.205]    0.10 0.180500
5    (0.205,0.254]    0.08 0.229500
6    (0.254,0.303]    0.03 0.278500

现在我们只需要绘制它。应该很容易...

ggplot(t_df, aes(means, density)) + 
  geom_histogram(stat = "identity") +
  geom_density(stat = "identity")

嗯，不是我想要的。可以肯定的是，我确实尝试过在 geom_density 中不使用 stat = "identity"，此时它抱怨没有 y。

#lets try adding ..density.. then
ggplot(t_df, aes(means, density)) + 
  geom_histogram(stat = "identity") +
  geom_density(aes(y = ..density..))

更奇怪。

好吧，也许让我们放弃从汇总数据中获取密度曲线。也许我们需要稍微混合一下这些方法......

#adding together
ggplot(t_df, aes(means, density)) +
  geom_bar(stat = "identity") +
  geom_density(data=t, aes(r, y = ..density..), stat = 'density')

好的，至少形状是现在的。现在，我们需要以某种方式缩小它。

#lets try dividing by the number of bins
ggplot(t_df, aes(means, density)) +
  geom_bar(stat = "identity") +
  geom_density(data=t, aes(r, y = ..density../20), stat = 'density')

看起来我们赢了。除了数字是硬编码的。

#removing the hardcoding?
divisor = nrow(t_df)
ggplot(t_df, aes(means, density)) +
  geom_bar(stat = "identity") +
  geom_density(data=t, aes(r, y = ..density../divisor), stat = 'density')

Error in eval(expr, envir, enclos) : object 'divisor' not found

嗯，我几乎预计它会起作用。现在我尝试在这里和那里添加一些..，还有..count..和..sum..，第一个给出了另一个错误的结果，第二个抛出了一个错误。我也尝试使用乘数（1/20），没有运气。

#salvation with get()
divisor = nrow(t_df)
ggplot(t_df, aes(means, density)) +
  geom_bar(stat = "identity") +
  geom_density(data=t, aes(r, y = ..density../get("divisor", pos = 1)), stat = 'density')

所以，我终于得到了正确的数字（我想；我希望）。

请告诉我有一种更简单的方法。

PS。 get() 技巧显然不适用于函数。我会在这里放置一个工作函数以供将来使用，但这也不是那么容易。

Answer 1

首先，阅读 Wickham 关于 R 中的密度，注意每个包/功能的弱点和特性。

密度总和为 1，但这并不意味着曲线/点不会超过 1。

与KernSmooth::bkde 相比，下面显示了这一点以及（至少）density 的默认值的不准确性（为简洁起见，使用基图）：

library(KernSmooth)
library(flux)
library(sfsmisc)

# uniform dist
set.seed(1)
dat <- runif(100)

d1 <- density(dat)
d1_ks <- bkde(dat)

par(mfrow=c(2,1))
plot(d1)
plot(d1_ks, type="l")

auc(d1$x, d1$y)
## [1] 1.000921

integrate.xy(d1$x, d1$y)
## [1] 1.000921

auc(d1_ks$x, d1_ks$y)
## [1] 1

integrate.xy(d1_ks$x, d1_ks$y)
## [1] 1

对 beta 版执行同样的操作：

# beta dist
set.seed(1)
dat <- rbeta(100, 0.5, 0.1)

d2 <- density(dat)
d2_ks <- bkde(dat)

par(mfrow=c(2,1))
plot(d2)
plot(d2_ks, typ="l")

auc(d2$x, d2$y)
## [1] 1.000187

integrate.xy(d2$x, d2$y)
## [1] 1.000188

auc(d2_ks$x, d2_ks$y)
## [1] 1

integrate.xy(d2_ks$x, d2_ks$y)
## [1] 1

auc 和 integrate.xy 都使用梯形规则，但我运行它们来显示这一点并显示两个不同函数的结果。

关键是密度实际上总和为 1，尽管 y 轴值使您相信它们不是。我不确定你想用你的操作解决什么问题。