I'm trying to compute the 95% HPDI of a posterior from 10000 draws from a distribution.
I've been instructed to use density() in R with default settings, and then sort the kernel density estimates in order to compute the HPDI.
From the documentation of the density() function:
The algorithm used in density.default disperses the mass of the empirical distribution function over a regular grid of at least 512 points and then uses the fast Fourier transform to convolve this approximation with a discretized version of the kernel and then uses linear approximation to evaluate the density at the specified points.
From the above section combined with the instruction I got I supposed it's the kernel$y values I should sort. Below is an illustration of my process
data = rnorm(10000) # Not the data I use, but works for the purpose of illustrating
kernel = density(data)
# Can plot the kernel and look at distribution
plot(kernel)
# From here a bit uncertain
sorted.mass = sort(kernel$y)
From here I don't get what to do. My initial thought was that I should now take the sum of sorted values untill I had 0.95, but sum(kernel$y) does not sum to one. Any input on how to proceed from here?
Edited according to Tims comment
data = rnorm(10000) # Not the data I use, but works for the purpose of illustrating
kernel = density(data)
# Can plot the kernel and look at distribution
plot(kernel)
# Numerically find HPDI
const = sum(kernel$y)
spxx = sort(kernel$y, decreasing = TRUE) / const
crit = spxx[which(cumsum(spxx) >=0.95)[1]]*const
abline(crit,0, col="blue")
