Jak znaleźć 95% wiarygodny przedział?

14

Próbuję obliczyć 95% wiarygodny przedział czasu następującego rozkładu tylnego. Nie mogłem znaleźć dla niej funkcji w R, ale czy poniższe podejście jest prawidłowe?

x <- seq(0.4,12,0.4)
px <-  c(0,0, 0, 0, 0, 0, 0.0002, 0.0037, 0.018, 0.06, 0.22 ,0.43, 0.64,0.7579, 0.7870, 0.72, 0.555, 0.37, 0.24, 0.11, 0.07, 0.02, 0.009, 0.005, 0.0001, 0,0.0002, 0, 0, 0)
plot(x,px, type="l")
mm <- sum(x*px)/sum(px)
var <- (sum((x)^2*px)/sum(px)) - (mm^2)
cat("95% credible interval: ", round(mm -1.96*sqrt(var),3), "-", round(mm + 1.96*sqrt(var),3),"\n")
użytkownik19758
źródło
1
Niezupełnie - przyjąłeś rozkład normalny i równy odstęp względem średniej, z których żadne nie jest szczególnie uzasadnione w tym kontekście. W rzeczywistości udało ci się uchwycić około prawdopodobieństwa, zakładając, że jest to rozkład dyskretny i musisz nieznacznie rozszerzyć interwał, aby uzyskać . Lepiej może być wziąć region o największej gęstości, który wynosi jeśli jest to rozkład dyskretny. Alternatywnie weź przedział czasu, aby prawdopodobieństwo, że będzie poniżej niego, wynosi lub mniej, a prawdopodobieństwo, że powyżej niego, wynosi lub mniej, również tutaj. 95 % [ 4,4 ; 8,0 ] 2,5 % 2,5 % [ 4,4 ; 8,0 ]94%95%[4.4,8.0]2.5%2.5%[4.4,8.0]
Henry

Odpowiedzi:

26

Jak zauważył Henry , zakładasz rozkład normalny i jest całkowicie w porządku, jeśli twoje dane są zgodne z rozkładem normalnym, ale będzie niepoprawny, jeśli nie możesz założyć dla niego rozkładu normalnego. Poniżej opisuję dwa różne podejścia, które można zastosować do nieznanego rozkładu, biorąc pod uwagę tylko punkty danych xi towarzyszące im szacunki gęstości px.

Pierwszą rzeczą do rozważenia jest to, co dokładnie chcesz podsumować za pomocą interwałów. Na przykład, możesz być zainteresowany interwałami uzyskanymi za pomocą kwantyli, ale możesz również być zainteresowany regionem o największej gęstości (zobacz tutaj lub tutaj ) swojej dystrybucji. Chociaż nie powinno to robić dużej różnicy (jeśli w ogóle) w prostych przypadkach, takich jak dystrybucje symetryczne, unimodalne, będzie to miało znaczenie dla bardziej „skomplikowanych” dystrybucji. Zasadniczo kwantyle podadzą przedział zawierający masę prawdopodobieństwa skoncentrowaną wokół mediany (środkowy twojego rozkładu), podczas gdy region o największej gęstości to obszar wokół trybów100α%dystrybucji. Będzie to wyraźniejsze, jeśli porównasz dwie wykresy na poniższym obrazku - kwantyle „wycinają” rozkład w pionie, a region o największej gęstości „wycina” go w poziomie.

Kwantyle vs przedziały HDR

Następną rzeczą do rozważenia jest sposób radzenia sobie z faktem, że masz niepełne informacje o rozkładzie (zakładając, że mówimy o ciągłym rozkładzie, masz tylko kilka punktów, a nie funkcję). Co możesz z tym zrobić, to wziąć wartości „takie, jakie są” lub użyć jakiegoś rodzaju interpolacji lub wygładzenia, aby uzyskać wartości „pomiędzy”.

Jednym podejściem byłoby użycie interpolacji liniowej (patrz ?approxfunR) lub alternatywnie coś bardziej gładkiego jak splajny (patrz ?splinefunR). Jeśli wybierzesz takie podejście, musisz pamiętać, że algorytmy interpolacji nie mają wiedzy domenowej o twoich danych i mogą zwracać nieprawidłowe wyniki, takie jak wartości poniżej zera itp.

# grid of points
xx <- seq(min(x), max(x), by = 0.001)

# interpolate function from the sample
fx <- splinefun(x, px) # interpolating function
pxx <- pmax(0, fx(xx)) # normalize so prob >0

Drugim podejściem, które można rozważyć, jest użycie rozkładu gęstości / mieszanki jądra w celu przybliżenia dystrybucji przy użyciu posiadanych danych. Trudność polega na tym, aby zdecydować o optymalnej przepustowości.

# density of kernel density/mixture distribution
dmix <- function(x, m, s, w) {
  k <- length(m)
  rowSums(vapply(1:k, function(j) w[j]*dnorm(x, m[j], s[j]), numeric(length(x))))
}

# approximate function using kernel density/mixture distribution
pxx <- dmix(xx, x, rep(0.4, length.out = length(x)), px) # bandwidth 0.4 chosen arbitrary

Następnie znajdziesz przedziały zainteresowania. Możesz postępować numerycznie lub symulacyjnie.

1a) Pobieranie próbek w celu uzyskania interwałów kwantylowych

# sample from the "empirical" distribution
samp <- sample(xx, 1e5, replace = TRUE, prob = pxx)

# or sample from kernel density
idx <- sample.int(length(x), 1e5, replace = TRUE, prob = px)
samp <- rnorm(1e5, x[idx], 0.4) # this is arbitrary sd

# and take sample quantiles
quantile(samp, c(0.05, 0.975)) 

1b) Pobieranie próbek w celu uzyskania regionu o największej gęstości

samp <- sample(pxx, 1e5, replace = TRUE, prob = pxx) # sample probabilities
crit <- quantile(samp, 0.05) # boundary for the lower 5% of probability mass

# values from the 95% highest density region
xx[pxx >= crit]

2a) Znajdź kwantyle numerycznie

cpxx <- cumsum(pxx) / sum(pxx)
xx[which(cpxx >= 0.025)[1]]   # lower boundary
xx[which(cpxx >= 0.975)[1]-1] # upper boundary

2b) Znajdź region o największej gęstości numerycznie

const <- sum(pxx)
spxx <- sort(pxx, decreasing = TRUE) / const
crit <- spxx[which(cumsum(spxx) >= 0.95)[1]] * const

Jak widać na poniższych wykresach, w przypadku unimodalnego, symetrycznego rozkładu obie metody zwracają ten sam interwał.

Dwa rodzaje interwałów

Oczywiście możesz także spróbować znaleźć interwał wokół jakiejś centralnej wartości, takiej jak i użyć pewnego rodzaju optymalizacji, aby znaleźć odpowiednią , ale dwa opisane powyżej podejścia wydają się być stosowane częściej i są bardziej intuicyjne.100α%Pr(Xμ±ζ)αζ

Tim
źródło
Dlaczego próbujesz, kiedy możesz po prostu obliczyć kwantyle bezpośrednio z podanych informacji (przy użyciu dowolnej metody)?
whuber
1
@ Whuber, ponieważ jest tani i łatwy, ale dokonam edycji, aby opisać jutro obliczenia inne niż symulacja.
Tim
Cześć Tim, to jest bardzo pomocne. Czy nie byłoby poprawne również wziąć kwantyl z rozkładu? (dolny <- x [który (logiczny (róż. (suma (px) / suma (px)> 0,025)))]) (górny <- x [który (logiczny (róż. (suma (px) / suma) (px) <0,975)))])
user19758
@ user19758 proszę sprawdzić moją edycję.
Tim
+1 Dodatkowe objaśnienia, ilustracje i kod ustanawiają wysoki standard odpowiedzi na tej stronie. Dziękuję Ci!
whuber