# How to find 95% credible interval?

I am trying to compute the 95% credible interval of the following posterior distribution. I could not find the function in R for it but is the approach below correct?

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

• Not really - you have assumed a normal distribution and an equal interval about the mean, neither of which are particularly justifiable in this context. In fact you have captured about $94\%$ of the probability, assuming this is a discrete distribution, and need to slightly widen your interval to get $95\%$. Better might be to take the highest density region which is $[4.4,8.0]$ if this is a discrete distribution. Alternatively take an interval so the probability of being below it is $2.5\%$ or less, and the probability of being above it is $2.5\%$ or less, also $[4.4,8.0]$ here. – Henry Oct 17 '16 at 20:30

As noted by Henry, you are assuming normal distribution and it's perfectly ok if your data follows normal distribution, but will be incorrect if you cannot assume normal distribution for it. Below I describe two different approaches that you could use for unknown distribution given only datapoints x and accompanying density estimates px.

The first thing to consider is what exactly do you want to summarize using your intervals. For example, you could be interested in intervals obtained using quantiles, but you could also be interested in highest density region (see here, or here) of your distribution. While this should not make much (if any) difference in simple cases like symmetric, unimodal distributions, this will make a difference for more "complicated" distributions. Generally, quantiles will give you interval containing probability mass concentrated around the median (the middle $100\alpha\%$ of your distribution), while highest density region is a region around the modes of the distribution. This will be more clear if you compare the two plots on the picture below -- quantiles "cut" the distribution vertically, while highest density region "cuts" it horizontally.

Next thing to consider is how to deal with the fact that you have incomplete information about the distribution (assuming that we are talking about continuous distribution, you have only a bunch of points rather then a function). What you could do about it is to take the values "as is", or use some kind of interpolation, or smoothing, to obtain the "in between" values.

One approach would be to use linear interpolation (see ?approxfun in R), or alternatively something more smooth like splines (see ?splinefun in R). If you choose such approach you have to remember that interpolation algorithms have no domain knowledge about your data and can return invalid results like values below zero 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


Second approach that you could consider is to use kernel density/mixture distribution to approximate your distribution using the data you have. The tricky part in here is to decide about optimal bandwidth.

# 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


Next, you are going to find the intervals of interest. You can either proceed numerically, or by simulation.

1a) Sampling to obtain quantile intervals

# 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) Sampling to obtain highest density region

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) Find quantiles numerically

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


2b) Find highest density region numerically

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


As you can see on the plots below, in case of unimodal, symmetric distribution both methods return the same interval.

Of course, you could also try to find $100\alpha\%$ interval around some central value such that $\Pr(X \in \mu \pm \zeta) \ge \alpha$ and use some kind of optimization to find appropriate $\zeta$, but the two approaches described above seem to be used more commonly and are more intuitive.

• Why do you sample when you could simply compute the quantiles directly from the information given (using either method)? – whuber Oct 17 '16 at 21:29
• @whuber because it's cheap and easy, but I'll edit to describe the non-simulation calculation tomorrow. – Tim Oct 17 '16 at 21:48
• Hi Tim, This is very helpful. Would it not be correct to also just take the quantile from the distibrution. (lower <- x[which(as.logical(diff(cumsum(px)/sum(px)>0.025)))]) (upper <- x[which(as.logical(diff(cumsum(px)/sum(px)<0.975)))]) – user19758 Oct 17 '16 at 23:44
• @user19758 please check my edit. – Tim Oct 18 '16 at 11:29
• +1 The additional explanations, illustrations, and code set a high standard for answers on this site. Thank you! – whuber Oct 18 '16 at 13:52