在两点之间着色核密度图.

JD *_*ong 92 plot r

我经常使用核密度图来说明分布.这些在R中创建简单快捷,如下所示:

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

这给了我这个漂亮的小PDF:

在此输入图像描述

我想将PDF下面的区域从第75百分位到第95百分位.使用quantile函数计算点很容易:

q75 <- quantile(draws, .75)
q95 <- quantile(draws, .95)
Run Code Online (Sandbox Code Playgroud)

但是我如何遮蔽q75和之间的区域q95

Dir*_*tel 73

有了这个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"))
Run Code Online (Sandbox Code Playgroud)

输出(由JDL添加)

在此输入图像描述

  • 如果你没有提供结构,我永远不会得到它.谢谢! (3认同)
  • 这是其中的一件事……自从黎明之前就一直存在于“演示(图形)”中,因此不时出现。NBER回归着色等的想法相同。 (2认同)

Ben*_*ker 69

另一种方案:

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)
Run Code Online (Sandbox Code Playgroud)

结果: 替代文字

  • 嘿那真棒!并充满了ggplot的善良! (2认同)

Mil*_*der 20

扩展的解决方案:

如果你想要遮蔽两个尾部(复制和粘贴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"))
Run Code Online (Sandbox Code Playgroud)

结果:

双尾聚


jor*_*ran 18

这个问题需要一个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))
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述