代码之家 › 专栏 › 技术社区 › JD Long

在两点之间绘制核密度图。

plot r

JD Long · 技术社区 · 15 年前

我经常用核密度图来说明分布。这些都是很容易和快速在R中创建的,就像这样:

set.seed(1)
绘制<-rnorm(100)^2
密度-密度(绘图)
情节(窝点)
#或者像这样的一行:绘图(密度(rnorm(100)^2))。
< /代码> 

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



我想把PDF下的区域从75%降到95%。使用分位数的功能很容易计算点:


q75<-分位数(draws,.75)
Q95<-分位数(绘图,.95)
< /代码> 

但如何在q75和q95之间对区域进行着色?


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



我想把PDF下的区域从75%降到95%。使用quantile功能:

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


但是我该怎么给中间的区域加阴影呢q75和q95?

5 回复 | 直到 7 年前

JD Long 8 年前

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


   
    
     您需要找到分位数值的索引,以获得实际值
     
      (x,y).
      
       pairs.
      
     
    
    
     
      
      
     
    
    
     
      
       
        编辑:
        
         开始吧:
        
       
      
     
    
    
     
      
       
        x1<-min(which(dens$x>=q75))。
x2<-最大值(dens$x<q95)
使用(dens,polygon(x=c(x[c(x1,x1:x2,x2)),y=c(0,y[x1:x2],0),col=“gray”))
< /代码> 

输出(由JDL添加)


分位数值的x,以获得实际值(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 15 年前

另一个解决方案:

dd<-with(dens,data.frame(x,y))。
图书馆(ggplot2)
qplot(x,y,data=dd,geom=“line”)。+
geom_功能区(数据=子集(dd,x>q75&x<q95)、aes(ymax=y)、ymin=0,
fill=“red”,colour=na,alpha=0.5)
< /代码> 

结果:

结果: alt text

zero323 little_kid_pea 10 年前

扩展的解决方案:

如果要对两个尾部进行着色(复制和粘贴dirk的代码),并使用已知的x值:

set.seed(1)
绘制<-rnorm(100)^2
密度-密度(绘图)
情节(窝点)

Q2和lt;- 2
问题65<-6.5
QN08<--0.8
问题02<--0.2

x1<-最小值(dens$x>=q2)
x2<-最大值(dens$x<q65)
x3<-最小值(dens$x>=qn08)
x4<-最大值(dens$x<qn02)

使用(dens,polygon(x=c(x[c(x1,x1:x2,x2)),y=c(0,y[x1:x2],0),col=“gray”))
使用(dens,polygon(x=c(x[c(x3,x3:x4,x4)),y=c(0,y[x3:x4],0),col=“gray”))
< /代码> 

结果:


ls(复制和粘贴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"))


结果:

joran 14 年前

这个问题需要一个 lattice answer。这是一个非常基本的方法,只需调整Dirk和其他人使用的方法:


   
    设置数据
SET种子(1)
绘制<-rnorm(100)^2
密度-密度(绘图)

#放入一个简单的数据框架
d<-data.frame(x=dens$x,y=dens$y)

#定义自定义面板功能;
#颜色等选项不需要硬编码
shadePanel<-函数(x,y,shadelims){
面板线条(x,y)
m1<-min(其中(x>=shadelims[1]))
m2<-最大值(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”)。
}

γ图
xyplot(y~x,data=d,panel=shadepanel,shadelims=c(1,3))。
< /代码>
    
     
    
    简单地采用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))

Matt Flor 7 年前

下面是另一个 ggplot2 variant,它基于一个函数,该函数近似于原始数据值下的内核密度:


   
    approvidens<-函数(x){
密度<-密度(x)
F<-带(dens,approxfun(x,y))。
f(x)
}
< /代码> 



使用原始数据(而不是使用密度估计的x和y值生成新的数据帧)还可以在分位数值取决于数据分组变量的分位数图中工作:



使用代码

library(tidyverse)
图书馆(RColorBrewer)

伪数据
SET种子(1)
N-LT;-1E2
dt<-tibble(值=rnorm(n)^2)

#在给定值处近似密度的函数
批准<-功能(X){
密度<-密度(x)
F<-带(dens,approxfun(x,y))。
f(x)
}

探针<-C(0.75,0.95)

dt<-dt%>%
突变(dy=approxdens(value),计算密度
P=百分比排名(值),百分比排名
pcat=as.因子(切割(P,断裂=探针,基于探针的百分比类别
include.lowest=true)))

ggplot(dt,aes(value,dy))。+
土工带(aes(ymin=0,ymax=dy,fill=pcat))。+
GeoMyLink()+
磅秤灌装机(guide=“none”)。+
主题BW()



#两组虚拟数据
DT2<-Tibble(类别=C(rep(“A”,N),rep(“B”,N)),
值=c(rnorm(n)^2,rnorm(n,mean=2)))

DT2<-DT2%>%
分组依据(类别)%>%
突变(dy=approxdens(value)
P=百分比排名(值)
pcat=as.因子(切割(P,断裂=探针,
include.lowest=true)))

刻面图
ggplot(dt2,aes(value,dy))。+
土工带(aes(ymin=0,ymax=dy,fill=pcat))。+
GeoMyLink()+
facet_wrap(~类别,nrow=2,scales=“fixed”)。+
磅秤灌装机(guide=“none”)。+
主题BW()
< /代码> 

由reprex package(v0.2.0)于2018-07-13创建。
t基于近似原始数据值的核密度的函数:

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()


创建于2018-07-13reprex package(v0.2.0)。