# Highest Posterior Density Interval (HPDI) using kernel

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