ggplot2直方图,密度曲线总和为1

Question

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

一些例子:

解决方案仅适用于标准化的正常数据 ggplot2:使用密度曲线叠加直方图

具有离散数据且无密度曲线 ggplot2密度直方图,宽度= .5,vline和居中条位置

没有答案 使用自定义分档使用ggplot2覆盖密度和直方图

在我的数据上,密度不总和为1 在ggplot2中创建密度直方图？

我的数据ggplot2密度直方图与自定义bin边缘不总和为1

这里用例子详细解释,但密度不是1,我的数据 "密度"曲线覆盖在直方图上,其中垂直轴是频率(即计数)或相对频率？

-

一些示例代码:

#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")

嗯,不是我想要的.可以肯定的是,我确实尝试过没有stat = "identity"geom_density,此时它抱怨没有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.

以下显示了这个和(至少)默认值的不准确性density,比如说,KernSmooth::bkde(为了简化输入,使用基础图):

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轴值导致您相信它们没有.我不确定你要用你的操作解决什么问题.