Comment trouver un intervalle crédible à 95%?

13

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")
user19758
la source
1
Pas vraiment - vous avez supposé une distribution normale et un intervalle égal autour de la moyenne, qui ne sont pas particulièrement justifiables dans ce contexte. En fait, vous avez capturé environ 94% de la probabilité, en supposant qu'il s'agit d'une distribution discrète, et vous devez légèrement élargir votre intervalle pour obtenir 95% . Mieux vaut peut-être prendre la région de densité la plus élevée qui est [4.4,8.0] s'il s'agit d'une distribution discrète. Vous pouvez également prendre un intervalle de sorte que la probabilité d'être en dessous soit de 2.5% ou moins, et la probabilité d'être au-dessus de lui est de 2.5% ou moins, également [4.4,8.0] ici.
Henry

Réponses:

24

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 modes100α%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.

Quantiles vs intervalles HDR

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.

Deux types d'intervalles

100α%ζPr(Xμ±ζ)αζ

Tim
la source
Pourquoi échantillonnez-vous alors que vous pouvez simplement calculer les quantiles directement à partir des informations fournies (en utilisant l'une ou l'autre méthode)?
whuber
1
@whuber parce que c'est bon marché et facile, mais je vais éditer pour décrire le calcul de non-simulation demain.
Tim
Salut Tim, c'est très utile. Ne serait-il pas juste de prendre également le quantile de la distribution. (inférieur <- x [qui (as.logical (diff (cumsum (px) / sum (px)> 0.025)))]) (supérieur <- x [which (as.logical (diff (cumsum (px) / sum) (px) <0,975)))])
user19758
@ user19758 veuillez vérifier ma modification.
Tim
+1 Les explications, illustrations et codes supplémentaires établissent un standard élevé pour les réponses sur ce site. Je vous remercie!
whuber