Feature request/idea: option to plot highest-density intervals

Now that arviz is an optional dependence, it would be nice to have the option to easily plot the highest-density intervals (from arviz.hdi) on the histograms rather than percentiles - in the case of asymmetric distributions, I find these more informative than the current quantiles behavior (though in the limit of a symmetric distribution they would be the same).

If you think this is a useful addition, I'd be happy to write up a PR. Before I do that, though, I thought I'd ask:

• Do you think this is useful?
• If so, I'd imagine it behaving similarly to quantiles by drawing lines on the 1-D histograms, though each probability value would give two lines, for the upper and lower bounds of that interval, rather than a single line.
• It could take an iterable of probabilities to draw, or just a single probability. If multiple probabilities were given, would it make sense to have different line styles by default?

Here's a quick example showing the difference with symmetric vs. asymmetric distributions - HDI 50% interval is in green, vs. 25% and 75% quantiles.

import corner
import numpy as np
import arviz as az 
import matplotlib.pyplot as plt

ndim, nsamples = 2, 20000
np.random.seed(42)
data =(np.random.randn(ndim * nsamples)).reshape([nsamples, ndim])

for samples in data, np.cos(data):
    figure = corner.corner(samples, quantiles=[0.25, 0.75])

    axes = np.array(figure.axes).reshape((ndim, ndim))

    # Loop over the diagonal
    for i in range(ndim):
        ax = axes[i, i]
        for prob in az.hdi(samples[:, i], hdi_prob=0.50):
            ax.axvline(prob, color="g")

Start at corner.corner and trace the existing quantiles behavior for one-dimensional histograms. Review how the optional arviz dependency is handled and clarify the open API questions in this issue; done means an option can use arviz HDI intervals and draw their bounds on the histograms.