R:找到密度图的最大值

Question

我有大约 25,000 行的数据 myData带列attr具有从 0 -> 45,600 的值。我不确定如何制作简化或可重复的数据...

无论如何，我正在绘制 attr 的密度像下面这样，我也找到了 attr密度最大的值:

library(ggplot)
max <- which.max(density(myData$attr)$y)
density(myData$attr)$x[max]
ggplot(myData, aes(x=attr))+ 
  geom_density(color="darkblue", fill="lightblue")+
  geom_vline(xintercept = density(myData$attr)$x[max])+
  xlab("attr")

这是我在最大点处得到的 x 截距图:


由于数据是倾斜的，然后我尝试通过添加 scale_x_log10() 以对数比例绘制 x 轴。到 ggplot ，这是新图:


我现在的问题是:

1. 为什么它现在有 2 个最高分？为什么我的 x 截距不再位于最大点？

2. 如何找到 2 个新的最大点的截距？

最后，我尝试将 y 轴转换为 count反而:

ggplot(myData, aes(x=attr)) +
  stat_density(aes(y=..count..), color="black", fill="blue", alpha=0.3)+
  xlab("attr")+
  scale_x_log10()

我得到了以下情节:


3. 我如何找到 count 2个峰？

Answer 1

为什么密度形状不同

为了将我的评论放在更完整的上下文中，ggplot 在进行密度估计之前获取日志，这会导致形状差异，因为分箱覆盖了域的不同部分。例如，

(bins <- seq(1, 10, length.out = 10))
#>  [1]  1  2  3  4  5  6  7  8  9 10
(bins_log <- 10^seq(log10(1), log10(10), length.out = 10))
#>  [1]  1.000000  1.291550  1.668101  2.154435  2.782559  3.593814  4.641589
#>  [8]  5.994843  7.742637 10.000000

library(ggplot2)

ggplot(data.frame(x = c(bins, bins_log), 
                  trans = rep(c('identity', 'log10'), each = 10)), 
       aes(x, y = trans, col = trans)) + 
    geom_point()

这种分箱会影响最终的密度形状。例如，比较一个未变换的密度:

d <- density(mtcars$disp)
plot(d)

到一个预先记录的:

d_log <- density(log10(mtcars$disp))
plot(d_log)

请注意，模式的高度会翻转!我相信您要的是第一个，但是在密度之后应用对数转换，即

d_x_log <- d
d_x_log$x <- log10(d_x_log$x)
plot(d_x_log)

这里的模式是相似的，只是被压缩了。

转移到 ggplot

转到 ggplot 时，要在对数转换之前进行密度估计，最容易在 ggplot 之外预先进行:

library(ggplot2)

d <- density(mtcars$disp)

ggplot(data.frame(x = d$x, y = d$y), aes(x, y)) + 
    geom_density(stat = "identity", fill = 'burlywood', alpha = 0.3) + 
    scale_x_log10()

寻找模式

当只有一个模式时，找到模式相对容易；只是 d$x[which.max(d$x)] .但是当您有多种模式时，这还不够好，因为它只会向您显示最高的一种。一个解决方案是有效地取导数并寻找斜率从正变为负的位置。我们可以用 diff 以数字方式做到这一点。 ，而且由于我们只关心结果是正还是负，调用 sign将所有内容都变成 -1 和 1。* 如果我们调用 diff在此基础上，除了最大值和最小值(分别为 -2 和 2)之外，所有内容都将为 0。然后我们可以查找 which值小于 0，我们可以使用它进行子集化。 (因为 diff 不会在末尾插入 NA s，您必须在索引中添加一个。)总而言之，旨在处理密度对象，

d <- density(mtcars$disp)

modes <- function(d){
    i <- which(diff(sign(diff(d$y))) < 0) + 1
    data.frame(x = d$x[i], y = d$y[i])
}

modes(d)
#>          x           y
#> 1 128.3295 0.003100294
#> 2 305.3759 0.002204658

d$x[which.max(d$y)]    # double-check
#> [1] 128.3295

我们可以将它们添加到我们的绘图中，它们会得到很好的转换:

ggplot(data.frame(x = d$x, y = d$y), aes(x, y)) + 
    geom_density(stat = "identity", fill = 'mistyrose', alpha = 0.3) + 
    geom_vline(xintercept = modes(d)$x) +
    scale_x_log10()

绘制计数而不是密度

要将 y 轴转换为计数而不是密度，请将 y 乘以观测数，该数存储在密度对象中为 n :

ggplot(data.frame(x = d$x, y = d$y * d$n), aes(x, y)) + 
    geom_density(stat = "identity", fill = 'thistle', alpha = 0.3) + 
    geom_vline(xintercept = modes(d)$x) +
    scale_x_log10()

在这种情况下，它看起来有点傻，因为只有 32 个观测值分布在一个广泛的域中，但是当 n 更大、域更小时，它更易于解释:

d <- density(diamonds$carat, n = 2048)

ggplot(data.frame(x = d$x, y = d$y * d$n), aes(x, y)) + 
    geom_density(stat = "identity", fill = 'papayawhip', alpha = 0.3) + 
    geom_point(data = modes(d), aes(y = y * d$n)) +
    scale_x_log10()

* 如果值正好是 0，则为 0，但这在这里不太可能，无论如何都可以正常工作。