Add multimodal HDI

[sethaxen]

This PR adds multimodal HDI and does not break any public functionality.

Internally, it introduces the following:

• a private interface for density estimation methods (only for HDI estimation), with implementations using KDE, histograms, and width 1 histograms for integer-valued data
• a private interface for HDI estimation methods, with implementations for our existing unimodal HDI as well as multimodal HDI using a density estimator
• a KDE implementation using reflection at data-derived bounds, which should probably be upstreamed to KernelDensity.jl
• an implementation of the Improved Sheather-Jones kernel bandwidth selector, which performs much better than Silverman's rule or the normal reference rule for multimodal distributions. This should probably also be upstreamed.
• a public API for selecting among HDI methods with specific presets

Fixes #7

Example

julia> using Distributions, PosteriorStats, Random

julia> Random.seed!(42);

julia> x = rand(MixtureModel([Normal(), Normal(300)]), 1_000);

julia> hdi(x; method=:multimodal)
2-element Vector{IntervalSets.ClosedInterval{Float64}}:
 -2.9324632517668725 .. 71.84082102773944
 227.51879889767156 .. 303.18936258853194

julia> hdi(x; method=:multimodal_sample)
2-element Vector{IntervalSets.ClosedInterval{Float64}}:
 -2.8852798738966734 .. 1.293563377013842
 296.6331884366964 .. 302.9035525499188

julia> x = rand(MixtureModel([Binomial(50, 0.2), Binomial(50, 0.8)]), 1_000);

julia> hdi(x; method=:multimodal)
2-element Vector{IntervalSets.ClosedInterval{Int64}}:
 5 .. 14
 35 .. 46

[sethaxen]

Supersedes #13