Comment trouver un intervalle crédible à 95%?

J'essaie de calculer l'intervalle crédible à 95% de la distribution postérieure suivante. Je n'ai pas pu trouver la fonction dans R pour cela mais l'approche ci-dessous est-elle correcte?

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

Comme l'a noté Henry , vous supposez une distribution normale et c'est parfaitement correct si vos données suivent une distribution normale, mais seront incorrectes si vous ne pouvez pas supposer une distribution normale pour elle. Ci-dessous, je décris deux approches différentes que vous pourriez utiliser pour une distribution inconnue étant donné uniquement les points de données xet les estimations de densité qui l'accompagnent px.

La première chose à considérer est ce que vous voulez résumer exactement en utilisant vos intervalles. Par exemple, vous pourriez être intéressé par les intervalles obtenus en utilisant des quantiles, mais vous pourriez également être intéressé par la région de densité la plus élevée (voir ici ou ici ) de votre distribution. Bien que cela ne devrait pas faire beaucoup de différence (le cas échéant) dans des cas simples comme les distributions symétriques et unimodales, cela fera une différence pour les distributions plus "compliquées". En général, les quantiles vous donneront un intervalle contenant une masse de probabilité concentrée autour de la médiane (les moyens de votre distribution), tandis que la région de densité la plus élevée est une région autour des modes$100\alpha\%$de la distribution. Cela sera plus clair si vous comparez les deux graphiques sur l'image ci-dessous - les quantiles "coupent" la distribution verticalement, tandis que la région de densité la plus élevée la "coupe" horizontalement.

La prochaine chose à considérer est de savoir comment traiter le fait que vous avez des informations incomplètes sur la distribution (en supposant que nous parlons de distribution continue, vous n'avez qu'un tas de points plutôt qu'une fonction). Ce que vous pourriez faire, c'est de prendre les valeurs "telles quelles", ou d'utiliser une sorte d'interpolation, ou de lissage, pour obtenir les valeurs "intermédiaires".

Une approche consisterait à utiliser une interpolation linéaire (voir ?approxfundans R), ou alternativement quelque chose de plus lisse comme des splines (voir ?splinefundans R). Si vous choisissez une telle approche, vous devez vous rappeler que les algorithmes d'interpolation n'ont aucune connaissance du domaine de vos données et peuvent renvoyer des résultats invalides comme des valeurs inférieures à zéro, 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

La deuxième approche que vous pourriez envisager est d'utiliser la distribution de densité / mélange du noyau pour approximer votre distribution en utilisant les données dont vous disposez. La partie délicate ici est de décider de la bande passante optimale.

# 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

Ensuite, vous allez trouver les intervalles d'intérêt. Vous pouvez soit procéder numériquement, soit par simulation.

1a) Échantillonnage pour obtenir des intervalles quantiles

# 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) Échantillonnage pour obtenir la région de densité la plus élevée

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) Trouver des quantiles numériquement

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

2b) Trouver numériquement la région de densité la plus élevée

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

Comme vous pouvez le voir sur les graphiques ci-dessous, en cas de distribution symétrique unimodale, les deux méthodes renvoient le même intervalle.

$100\alpha\%$$\Pr(X \in \mu \pm \zeta) \ge \alpha$$\zeta$