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?
As pointed out in the other answers, 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".
This method has some advantages (and is preferred by many Bayesian statisticians for summarizing posterior densities), but unlike the quantile-based approaches is not invariant to transformations of the variable.
library(coda)
HPDinterval(as.mcmc(x), prob=0.9)
R implementation giving a $O((1-p)n + n\log(n))$ algorithm for probability $p$ and data size $n:$ k <- ceiling(0.9 * length(x)) - 1; x[c(0, k) + which.min(diff(sort(x), k))]
$\endgroup$
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
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.
