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

9

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?

ECII
źródło
Twoje pytanie jest trochę niejasne na temat „przedziałów, w których ...” - może być wiele przedziałów. Czy jesteś zainteresowany wyłącznie wewnętrznym 90%, tj. Symetrycznym przycinaniem z każdej strony? W końcu, od minimum do 90% ile, 90% danych jest przechwytywanych, podobnie przez 10% ile do wartości maksymalnej.
Iterator
Czy szukasz najkrótszego przedziału, przedziału symetrycznego (równe prawdopodobieństwo na każdym końcu), czy czegoś innego?
Glen_b

Odpowiedzi:

19

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)
Ben Bolker
źródło
3

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.

DWin
źródło
3

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
James
źródło
-1

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

Carl Witthoft
źródło