undesirable behaviour of plot_hdi, and a proposal

[drbenvincent]

Ok, so I have noticed behaviour of az.plot_hdi which I do not believe is a bug as such, but does represent undesirable behaviour.

A frequent use of arviz is to plot hdi intervals of a function for which you have a set of MCMC samples. Consider this example...

import arviz as az
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

RANDOM_SEED = 1234
rng = np.random.default_rng(RANDOM_SEED)

def f(x, s, c):
    return s * np.tanh(x/(s*c))
    
s = xr.DataArray(data=rng.beta(a=1, b=20, size=(4, 1000)), dims=['chain', 'draw'])
c = xr.DataArray(data=rng.normal(loc=8, scale=3, size=(4, 1000)), dims=['chain', 'draw'])
x = xr.DataArray(data=np.linspace(0, 1, 500), dims='x')

Y = f(x, s, c)

fig, ax =plt.subplots()
az.plot_hdi(x, Y, hdi_prob=0.95, smooth=False, ax=ax)
az.plot_hdi(x, Y, hdi_prob=0.50, smooth=False, ax=ax)
ax.plot(x, Y.median(dim=("chain", "draw")));

Which results in this.

There is a certain point in the x-axis where the hdi interval shifts (around x=0.5). I believe this is because the underlying distribution is asymmetric, and so you can get these discontinuities. So while this is not a bug (it does represent the actual hdi), it doesn't represent desirable behaviour.

So how about a new plot function which derives from plot_hdi, but which plots percentiles?

[twiecki]

I think Michael Betancourt also argued against HDIs because they require an estimator over the posterior.

[michaelosthege]

Maybe try pm.gp.util.plot_gp_dist? It creates something looking like a density based on percentiles, but I'd expect that it also creates a misleading result.

Instead, drawing a pm.gp.util.plot_gp_dist-like density band but based on KDEs might be more accurate.

You can always overlay HDI bounds or individual percentiles as a dashed line if needed.

[lucianopaz]

We could also use quantiles (I thought that arviz did that by default):

import arviz as az
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

RANDOM_SEED = 1234
rng = np.random.default_rng(RANDOM_SEED)

def f(x, s, c):
    return s * np.tanh(x/(s*c))
    
s = xr.DataArray(data=rng.beta(a=1, b=20, size=(4, 1000)), dims=['chain', 'draw'])
c = xr.DataArray(data=rng.normal(loc=8, scale=3, size=(4, 1000)), dims=['chain', 'draw'])
x = xr.DataArray(data=np.linspace(0, 1, 500), dims='x')

Y = f(x, s, c)

hdi_50 = Y.quantile([0.25, 0.75], dim=["chain", "draw"])
hdi_95 = Y.quantile([0.025, 0.975], dim=["chain", "draw"])

fig, ax =plt.subplots()
ax.fill_between(x, hdi_95.isel(quantile=0), hdi_95.isel(quantile=1), color="C1", alpha=0.5)
ax.fill_between(x, hdi_50.isel(quantile=0), hdi_50.isel(quantile=1), color="C1", alpha=0.5)
ax.plot(x, Y.median(dim=("chain", "draw")))

[ahartikainen]

And if sample size is 100?

[lucianopaz]

With 100 sample size I get this plot using the quantiles:

And this plot using plot_hdi:

[cluhmann]

This is the sort of argument that people use to motivate equal-tailed intervals, which are basically just the quantiles that @lucianopaz described. If you poke around, you can find other odd, possibly-undesirable behavior with HDIs (e.g., multi-modality) even though the HDI is always "correct".

[OriolAbril]

As an aside, az.hdi has a multimodal argument to find "discontinuous" hdi intervals around the multiple modes. It also has a circular argument (case in which using medians or equal-tailed intervals generally makes no sense at all)

[OriolAbril]

You can already use other estimators for the probability ranges, even if the name seems to indicate the opposite. I do agree we should make it easier or maybe even default to equal tail intervals, but there will always be arguments for both in my opinion and the choice should probably depend on the problem/model.

We currently have #1998 open to use the median and eti in the summary for example.

As example, change the plotting code to:

quantiles = Y.quantile((.025, .25, .5, .75, .975), dim=("chain", "draw")).transpose()

fig, ax =plt.subplots()
az.plot_hdi(x, hdi_data=quantiles.sel(quantile=[.025, .975]), ax=ax, fill_kwargs={"label": "95% ETI"})
az.plot_hdi(x, hdi_data=quantiles.sel(quantile=[.25, .75]), ax=ax, fill_kwargs={"label": "50% ETI"})
ax.plot(x, quantiles.sel(quantile=.5), label="Median")
ax.legend();

so equal tail intervals are used via hdi so you can even take advantage of the smooting if you want. And you can always use fill_between directly (it is what plot_hdi uses) also with any interval estimator function.

[drbenvincent]

Thanks for everyone's input on this. I think it could be worth considering @OriolAbril 's suggestion of changing default behaviour. Maybe plot_hdi could be renamed (in a slow deprication phase out) to something like plot_intervals?