A simple package for computing the highest density region of a univariate distribution defined in Distributions.jl. It is only intended for use on unimodal distributions as the package assumes that there is a single, connected, highest density region. Both continuous and discrete distributions are supported. However, the assumptions of the method may break down for discrete distributions; the method of O'Neill (2022) may be more appropriate. This can be seen in the occasional inconsistency by width 1 between the HDR found by the grid-based and optimization-based approaches. The exported function will not error on bimodal distributions, but it will not identify the correct highest density regions.

A grid-based approach and optimization-based approach are implemented. The optimisation approach will, in general, require fewer distribution quantile evaluations for the same level of accuracy. However, it requires loading Optimization.jl and therefore requires more memory.

This is a performant alternative to HighestDensityRegions.jl when the distribution of interest is univariate and unimodal.

To install the package, use the following command inside the Julia REPL:

using Pkg
Pkg.add("UnivariateUnimodalHighestDensityRegion")

To load the package, use the command:

using UnivariateUnimodalHighestDensityRegion

There is a single exported function, univariate_unimodal_HDR, which is used with univariate distributions from Distributions.jl.

After loading the package using the previous command we can find the highest density region of univariate distributions. Finding the 95% HDR of a Normal distribution will return the 2.5% and 97.5% quantiles; note, the distribution is symmetric so the method is unnecessary.

univariate_unimodal_HDR(Normal(0,2), 0.95)

2-element MVector{2, Float64} with indices SOneTo(2):
 -3.919927969080115
  3.919927969080115

The function is most valuable for asymmetric distributions such as a LogNormal or Poisson distribution:

univariate_unimodal_HDR(LogNormal(1,0.5), 0.95)

2-element MVector{2, Float64} with indices SOneTo(2):
 0.721779994018427
 6.312112357076725

univariate_unimodal_HDR(Poisson(4), 0.95)

2-element MVector{2, Int64} with indices SOneTo(2):
 1
 8

Repeating the previous examples with the optimization approach requires loading Optimization.jl and an additional package that contains optimization algorithms. Here we use OptimizationNLopt.jl. If they have not yet been installed they will also need to be installed using Pkg.add.

using Pkg
Pkg.add(["Optimization", "OptimizationNLopt"])
using Optimization, OptimizationNLopt
solver = NLopt.LN_BOBYQA()

Normal:

univariate_unimodal_HDR(Normal(0,2), 0.95, solver)

2-element MVector{2, Float64} with indices SOneTo(2):
 -3.919927969080115
  3.919927969080115

LogNormal:

univariate_unimodal_HDR(LogNormal(1,0.5), 0.95, solver)

2-element MVector{2, Float64} with indices SOneTo(2):
 0.7104607462933641
 6.300409489807155

Poisson:

univariate_unimodal_HDR(Poisson(4), 0.95, solver)

2-element MVector{2, Int64} with indices SOneTo(2):
 1
 8