Znajdź przedziały gęstości prawdopodobieństwa

Mam wektor

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)

(mój rzeczywisty wektor ma długość> 10 000) i chciałbym znaleźć przedziały, w których leży 90% gęstości. Czy quantile(x, probs=c(0.05,0.95), type=5)jest najbardziej odpowiedni, czy jest jakiś inny sposób?

Jak wskazano powyżej, istnieje wiele różnych sposobów definiowania przedziału obejmującego 90% gęstości. Tym, który nie został jeszcze wskazany, jest najwyższy [tylny] przedział gęstości ( wikipedia ), który jest zdefiniowany jako „najkrótszy przedział, dla którego różnica w wartościach empirycznej funkcji gęstości skumulowanej punktów końcowych jest nominalnym prawdopodobieństwem”.

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

Z pewnością wydaje się to najprostszym podejściem. Funkcja jest dość szybka. Używam go cały czas na próbkach, które są setki razy większe niż ta, której używasz, a stabilność szacunków powinna być dobra przy wielkości twojej próbki.

Istnieją funkcje w innych pakietach, które zapewniają bardziej kompletne zestawy statystyk opisowych. Ten, którego używam, jest Hmisc::describe, ale istnieje kilka innych pakietów z describefunkcjami.

Twoja droga wydaje się sensowna, szczególnie w przypadku dyskretnych danych w przykładzie,

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

ale innym sposobem byłoby użycie jądra o gęstości obliczeniowej:

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

Tak. :-). Być może wynik stats::densitybędzie bardziej użyteczny.