¿Cómo encontrar el intervalo creíble del 95%?

Estoy tratando de calcular el intervalo creíble del 95% de la siguiente distribución posterior. No pude encontrar la función en R para ello, pero ¿es correcto el siguiente enfoque?

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

Como señaló Henry , está asumiendo una distribución normal y está perfectamente bien si sus datos siguen una distribución normal, pero será incorrecta si no puede asumir una distribución normal para ella. A continuación, describo dos enfoques diferentes que podría usar para una distribución desconocida dados solo los puntos de datos xy las estimaciones de densidad que lo acompañan px.

Lo primero que debe considerar es qué es exactamente lo que desea resumir con sus intervalos. Por ejemplo, podría estar interesado en los intervalos obtenidos utilizando cuantiles, pero también podría estar interesado en la región de mayor densidad (ver aquí o aquí ) de su distribución. Si bien esto no debería hacer mucha (si alguna) diferencia en casos simples como distribuciones simétricas y unimodales, esto hará una diferencia para distribuciones más "complicadas". En general, los cuantiles le proporcionarán un intervalo que contiene la masa de probabilidad concentrada alrededor de la mediana (el de su distribución), mientras que la región de mayor densidad es una región alrededor de los modos$100\alpha\%$de la distribución. Esto será más claro si compara las dos gráficas en la imagen a continuación: los cuantiles "cortan" la distribución verticalmente, mientras que la región de mayor densidad la "corta" horizontalmente.

Lo siguiente a considerar es cómo lidiar con el hecho de que tiene información incompleta sobre la distribución (suponiendo que estamos hablando de distribución continua, solo tiene un montón de puntos en lugar de una función). Lo que podría hacer al respecto es tomar los valores "tal cual", o utilizar algún tipo de interpolación, o suavizado, para obtener los valores "intermedios".

Un enfoque sería usar interpolación lineal (ver ?approxfunen R), o alternativamente algo más suave como las splines (ver ?splinefunen R). Si elige este enfoque, debe recordar que los algoritmos de interpolación no tienen conocimiento de dominio sobre sus datos y pueden devolver resultados no válidos, como valores por debajo de cero, etc.

# 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

El segundo enfoque que podría considerar es utilizar la distribución de la densidad / mezcla del núcleo para aproximar su distribución utilizando los datos que tiene. La parte difícil aquí es decidir sobre el ancho de banda óptimo.

# 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

A continuación, encontrará los intervalos de interés. Puede proceder numéricamente o por simulación.

1a) Muestreo para obtener intervalos cuantiles

# 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) Muestreo para obtener la región de mayor densidad

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) Encuentra cuantiles numéricamente

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

2b) Encuentra la región de mayor densidad numéricamente

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

Como puede ver en los gráficos a continuación, en caso de distribución simétrica unimodal, ambos métodos devuelven el mismo intervalo.

Por supuesto, también podría intentar encontrar el intervalo alrededor de algún valor central tal que y usar algún tipo de optimización para encontrar apropiado , pero los dos enfoques descritos anteriormente parecen usarse más comúnmente y son más intuitivos.$100\alpha\%$$\Pr(X \in \mu \pm \zeta) \ge \alpha$$\zeta$