views:

233

answers:

1

I frequently use kernel density plots to illustrate distributions. These are easy and fast to create in R like so:

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

Which gives me this nice little PDF:

It's my PDF, not Adobe's

I'd like to shade the area under the PDF from the 75th to 95th percentiles. It's easy to calculate the points using the quantile function:

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

But how do I shade the the area between q75 and q95?

+10  A: 

With the polygon() function, see its help page and I believe we had similar questions here too.

You need to find the index of the quantile values to get the actual (x,y) pairs.

Edit: Here you go:

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

Output (added by JDL)

alt text

Dirk Eddelbuettel
I never would have gotten that to work if you had not provided the structure. Thanks!
JD Long
It's one of those things ... that have been in `demo(graphics)` since before the dawn on time so one comes across every now and then. Same idea for NBER regression shading etc.
Dirk Eddelbuettel
ohhhh. I KNEW I had seen it somewhere but could not pull from my mental index where I had seen it. I'm glad your mental index is better than mine.
JD Long
Thanks for the updated chart!
Dirk Eddelbuettel