1
$\begingroup$

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

enter image description here

$\endgroup$

1 Answer 1

1
$\begingroup$

The probability densities are not probabilities, don't sum to 1, they integrate to 1. For how to find the interval, you can check the How to find 95% credible interval? thread that includes the code.

$\endgroup$
3
  • $\begingroup$ Hi Tim, thanks for your response. It might be the plot that confuses me a bit but I don't understand how you use the critical value when finding highest density numerically. I'm able to recreate your process, and added an abline to illustrate it. But HPDI is intendend to give us a range between two (or more depending on shape) values which marks the HPDI. In your plot you seem to provide the boundary but not the values. It can be that I don't quite understand the process of what you're doing, but would you care to elaborate? $\endgroup$
    – OLGJ
    Commented Apr 7, 2022 at 15:57
  • $\begingroup$ @OLGJ as said in the linked answer, if the distribution is multimodal it doesn't have to be range. Unless you mean quantile interval, as mentioned. $\endgroup$
    – Tim
    Commented Apr 7, 2022 at 16:03
  • $\begingroup$ I meant that in the plot you linked to How to find 95% credible interval where HPDI is plotted, the abline used separates the distribution in a top/bottom region, and the values between where the line is tangent to the distribution marks the HPDI area. So in that sense the plot is confusing to me as it does not contain the HPDI area, but rather it highlights a separation line? $\endgroup$
    – OLGJ
    Commented Apr 8, 2022 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.