在两点之间对内核密度图进行阴影处理。

我经常使用内核密度图来说明分布。可以像这样在R中轻松便捷地创建它们：

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)
#or in one line like this: plot(density(rnorm(100)^2))

这给了我这个小巧的PDF：

我想将PDF下的区域的阴影从第75个百分位数增加到第95个百分位数。使用quantile函数很容易计算点：

q75 <- quantile(draws, .75)
q95 <- quantile(draws, .95)

但是如何给q75和q95之间的区域加阴影呢？

您能否提供在范围外部与范围内部之间加阴影的示例？谢谢。

#1 楼

使用polygon()函数，请参见其帮助页面，我相信这里也有类似的问题。

您需要找到分位数的索引才能获取实际的(x,y)对。

编辑：在这里，您可以进行以下操作：

x1 <- min(which(dens$x >= q75))  
x2 <- max(which(dens$x <  q95))
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))

输出（由JDL添加）

如果您没有提供该结构，我永远不会使它起作用。谢谢！

–JD Long
10年8月16日在17:17

这是其中的一件事……自从黎明前准时出现在演示（图形）中以来，时不时出现。 NBER回归着色等的想法相同。

–德克·埃德尔布特尔（Dirk Eddelbuettel）
10年8月16日在17:19

哦。我知道我在某处看到过它，但无法从我看到它的心理指数中得出。我很高兴您的心理指数比我的更好。

–JD Long
10年8月16日在17:20

#2 楼

另一个解决方案：

dd <- with(dens,data.frame(x,y))

library(ggplot2)

qplot(x,y,data=dd,geom="line")+
  geom_ribbon(data=subset(dd,x>q75 & x<q95),aes(ymax=y),ymin=0,
              fill="red",colour=NA,alpha=0.5)

结果：

#3 楼

扩展的解决方案：

如果要阴影化两条尾巴（复制并粘贴Dirk的代码）并使用已知的x值：

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)

q2     <- 2
q65    <- 6.5
qn08   <- -0.8
qn02   <- -0.2

x1 <- min(which(dens$x >= q2))  
x2 <- max(which(dens$x <  q65))
x3 <- min(which(dens$x >= qn08))  
x4 <- max(which(dens$x <  qn02))

with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))
with(dens, polygon(x=c(x[c(x3,x3:x4,x4)]), y= c(0, y[x3:x4], 0), col="gray"))

结果：

我有png文件并将其托管在freeimagehosting上，并且可能无法加载，因为...我不确定。

–Milktrader
2011-3-25在17:55

文件非常模糊。您能否重新创建并直接将其上传到此处，所以为此拥有自己的服务器服务？

–德克·埃德尔布特尔（Dirk Eddelbuettel）
2011年3月26日18:27

很抱歉，但是我看不到如何直接将其上传到SO。

–Milktrader
2011-3-28的1:03

#4 楼

这个问题需要一个lattice答案。这是一个非常基本的方法，只需修改Dirk等人使用的方法即可：

#Set up the data
set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)

#Put in a simple data frame   
d <- data.frame(x = dens$x, y = dens$y)

#Define a custom panel function;
# Options like color don't need to be hard coded    
shadePanel <- function(x,y,shadeLims){
    panel.lines(x,y)
    m1 <- min(which(x >= shadeLims[1]))
    m2 <- max(which(x <= shadeLims[2]))
    tmp <- data.frame(x1 = x[c(m1,m1:m2,m2)], y1 = c(0,y[m1:m2],0))
    panel.polygon(tmp$x1,tmp$y1,col = "blue")
}

#Plot
xyplot(y~x,data = d, panel = shadePanel, shadeLims = c(1,3))

#5 楼

这是基于函数的另一个ggplot2变体，该函数在原始数据值处近似内核密度：

 approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

使用原始数据（而不是使用密度估计的x和y值生成新的数据框）还可以在分位数值取决于变量的多面图中工作分组数据的依据：

使用的代码

 library(tidyverse)
library(RColorBrewer)

# dummy data
set.seed(1)
n <- 1e2
dt <- tibble(value = rnorm(n)^2)

# function that approximates the density at the provided values
approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

probs <- c(0.75, 0.95)

dt <- dt %>%
    mutate(dy = approxdens(value),                         # calculate density
           p = percent_rank(value),                        # percentile rank 
           pcat = as.factor(cut(p, breaks = probs,         # percentile category based on probs
                                include.lowest = TRUE)))

ggplot(dt, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    scale_fill_brewer(guide = "none") +
    theme_bw()



# dummy data with 2 groups
dt2 <- tibble(category = c(rep("A", n), rep("B", n)),
              value = c(rnorm(n)^2, rnorm(n, mean = 2)))

dt2 <- dt2 %>%
    group_by(category) %>% 
    mutate(dy = approxdens(value),    
           p = percent_rank(value),
           pcat = as.factor(cut(p, breaks = probs,
                                include.lowest = TRUE)))

# faceted plot
ggplot(dt2, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    facet_wrap(~ category, nrow = 2, scales = "fixed") +
    scale_fill_brewer(guide = "none") +
    theme_bw()

由reprex包（v0.2.0）于2018-07-13创建。

编程黑洞网