Jak znaleźć 95% wiarygodny przedział?

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")

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ów$100\alpha\%$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.

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ł.

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\alpha\%$$\Pr(X \in \mu \pm \zeta) \ge \alpha$$\zeta$