# Find probability density intervals

I have the vector

x <- c(1,2,3,4,5,5,5,6,6,6,6,
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
7,7,7,7,7,7,7,7,8,8,8,8,9,9,9,10)


(my actual vector has a length of >10,000) and I would like to find the intervals where the 90% of the density lies. Is quantile(x, probs=c(0.05,0.95), type=5) the most appropriate or is there any other way?

## migrated from stackoverflow.comNov 17 '11 at 16:31

This question came from our site for professional and enthusiast programmers.

• Your question is a little vague about "the intervals where..." - there could be multiple intervals. Are you interested in solely the inner 90%, i.e. symmetrically trimming on each side? After all, from the minimum to 90%ile, 90% of the data is captured, similarly for 10%ile to the max value. – Iterator Nov 16 '11 at 12:53
• Are you looking for a shortest interval, a symmetric interval (equal probability out each end), or something else? – Glen_b Sep 3 '16 at 3:14

As pointed out above, there are many different ways to define an interval that includes 90% of the density. One that hasn't been pointed out yet is the highest [posterior] density interval (wikipedia), which is defined as "the shortest interval for which the difference in the empirical cumulative density function values of the endpoints is the nominal probability".

library(coda)
HPDinterval(as.mcmc(x), prob=0.9)


It certainly seems like the most straightforward approach. The function is quite fast. I use it all the time on samples that are hundreds of times larger that the one you are using, and the stability of the estimates should be good at your sample size.

There are functions in other packages that provide more complete sets of descriptive statistics. The one I use is Hmisc::describe, but there are several other packages with describe functions.

Your way seems sensible, especially with the discrete data in the example,

quantile(x,probs=c(0.05,0.95), type=5)
5% 95%
2.8 9.0


but another way would be to use a computed density kernel:

dx <- density(x)
dn <- cumsum(dx$y)/sum(dx$y)
li <- which(dn>=0.05)[1]
ui <- which(dn>=0.95)[1]
dx\$x[c(li,ui)]
[1] 2.787912 9.163246


Yes. :-). You may find the output of stats::density to be more useful.