Sombreando un gráfico de densidad de kernel entre dos puntos.

Frecuentemente uso diagtwigs de densidad del kernel para ilustrar las distribuciones. Estos son fáciles y rápidos de crear en R así que:

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

Lo que me da este pequeño y bonito PDF:

enter image description here

Me gustaría sombrear el área debajo del PDF del percentil 75 al 95. Es fácil calcular los puntos usando la función de quantile :

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

¿Pero cómo sombreo el área entre q75 y q95 ?

Con la función de polygon() , consulte su página de ayuda y creo que también tuvimos preguntas similares aquí.

Necesita encontrar el índice de los valores cuantílicos para obtener los pares reales (x,y) .

Editar: Aquí tienes:

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

Salida (agregada por JDL)

enter image description here

Otra solución:

 dd <- with(dens,data.frame(x,y)) library(ggplot2) qplot(x,y,data=dd,geom="line")+ geom_ribbon(data=subset(dd,x>q75 & x 

Resultado: texto alternativo

Una solución expandida:

Si quería sombrear ambas colas (copiar y pegar del código de Dirk) y usar valores x conocidos:

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

Resultado:

Polietileno de 2 colas

Esta pregunta necesita una respuesta lattice . Aquí hay uno muy básico, simplemente adaptando el método empleado por Dirk y otros:

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

enter image description here

Aquí hay otra variante de ggplot2 basada en una función que aproxima la densidad del kernel en los valores de datos originales:

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

Usar los datos originales (en lugar de producir un nuevo dataframe con los valores x e y de la estimación de densidad) tiene el beneficio de trabajar también en gráficos facetados donde los valores de los cuantiles dependen de la variable por la cual se agrupan los datos:

Código utilizado

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

Creado en 2018-07-13 por el paquete reprex (v0.2.0).