Implement multimodal HDI

[sethaxen]

Our current HDI estimator hdi assumes the distribution is unimodal. To support multimodal distributions, we would first fit a KDE (or histogram for discrete points) to the draws and then find the HDI from that. This is mostly useful for plotting.

Once #6 is complete, this could probably be easily done with https://github.com/tpapp/HighestDensityRegions.jl. See also Python ArviZ's implementation.

[sethaxen]

We should also look into estimating the HDI from average shifted histograms (see e.g. https://github.com/joshday/AverageShiftedHistograms.jl), as these are much faster to compute than the KDE. Probably a simulation study is necessary to compare both performance and accuracy of estimating known HDIs from KDE vs ASH.