"""Plot highest density intervals for regression data."""
import warnings
import numpy as np
from scipy.interpolate import griddata
from scipy.signal import savgol_filter
from xarray import Dataset
from ..rcparams import rcParams
from ..stats import hdi
from .plot_utils import get_plotting_function
[docs]
def plot_hdi(
x,
y=None,
hdi_prob=None,
hdi_data=None,
color="C1",
circular=False,
smooth=True,
smooth_kwargs=None,
figsize=None,
fill_kwargs=None,
plot_kwargs=None,
hdi_kwargs=None,
ax=None,
backend=None,
backend_kwargs=None,
show=None,
):
r"""Plot HDI intervals for regression data.
Parameters
----------
x : array-like
Values to plot.
y : array-like, optional
Values from which to compute the HDI. Assumed shape ``(chain, draw, \*shape)``.
Only optional if ``hdi_data`` is present.
hdi_data : array_like, optional
Precomputed HDI values to use. Assumed shape is ``(*x.shape, 2)``.
hdi_prob : float, optional
Probability for the highest density interval. Defaults to ``stats.hdi_prob`` rcParam.
See :ref:`this section <common_ hdi_prob>` for usage examples.
color : str, default "C1"
Color used for the limits of the HDI and fill. Should be a valid matplotlib color.
circular : bool, default False
Whether to compute the HDI taking into account ``x`` is a circular variable
(in the range [-np.pi, np.pi]) or not. Defaults to False (i.e non-circular variables).
smooth : boolean, default True
If True the result will be smoothed by first computing a linear interpolation of the data
over a regular grid and then applying the Savitzky-Golay filter to the interpolated data.
smooth_kwargs : dict, optional
Additional keywords modifying the Savitzky-Golay filter. See
:func:`scipy:scipy.signal.savgol_filter` for details.
figsize : (float, float), optional
Figure size. If ``None``, it will be defined automatically.
fill_kwargs : dict, optional
Keywords passed to :meth:`mpl:matplotlib.axes.Axes.fill_between`
(use ``fill_kwargs={'alpha': 0}`` to disable fill) or to
:meth:`bokeh.plotting.Figure.patch`.
plot_kwargs : dict, optional
HDI limits keyword arguments, passed to :meth:`mpl:matplotlib.axes.Axes.plot` or
:meth:`bokeh.plotting.Figure.patch`.
hdi_kwargs : dict, optional
Keyword arguments passed to :func:`~arviz.hdi`. Ignored if ``hdi_data`` is present.
ax : axes, optional
Matplotlib axes or bokeh figures.
backend : {"matplotlib", "bokeh"}, default "matplotlib"
Select plotting backend.
backend_kwargs : dict, optional
These are kwargs specific to the backend being used, passed to
:func:`matplotlib.pyplot.subplots` or :class:`bokeh.plotting.figure`.
For additional documentation check the plotting method of the backend.
show : bool, optional
Call backend show function.
Returns
-------
axes : matplotlib axes or bokeh figures
See Also
--------
hdi : Calculate highest density interval (HDI) of array for given probability.
Examples
--------
Plot HDI interval of simulated random-walk data using `y` argument:
.. plot::
:context: close-figs
>>> import numpy as np
>>> import arviz as az
>>> # time-steps random walk
>>> x_data =np.arange(0,100)
>>> # Mean random walk
>>> mu = np.zeros(100)
>>> for i in x_data: mu[i] = mu[i-1] + np.random.normal(0, 1, 1)
>>> # Simulated pp samples form the random walk time series
>>> y_data = np.random.normal(2 + mu * 0.5, 0.5, size = (2, 50, 100))
>>> az.plot_hdi(x_data, y_data)
``plot_hdi`` can also be given precalculated values with the argument ``hdi_data``. This example
shows how to use :func:`~arviz.hdi` to precalculate the values and pass these values to
``plot_hdi``. Similarly to an example in ``hdi`` we are using the ``input_core_dims``
argument of :func:`~arviz.wrap_xarray_ufunc` to manually define the dimensions over which
to calculate the HDI.
.. plot::
:context: close-figs
>>> hdi_data = az.hdi(y_data, input_core_dims=[["draw"]])
>>> ax = az.plot_hdi(x_data, hdi_data=hdi_data[0], color="r", fill_kwargs={"alpha": .2})
>>> az.plot_hdi(x_data, hdi_data=hdi_data[1], color="k", ax=ax, fill_kwargs={"alpha": .2})
``plot_hdi`` can also be used with Inference Data objects. Here we use the posterior predictive
to plot the HDI interval.
.. plot::
:context: close-figs
>>> X = np.random.normal(0,1,100)
>>> Y = np.random.normal(2 + X * 0.5, 0.5, size=(2,10,100))
>>> idata = az.from_dict(posterior={"y": Y}, constant_data={"x":X})
>>> x_data = idata.constant_data.x
>>> y_data = idata.posterior.y
>>> az.plot_hdi(x_data, y_data)
"""
if hdi_kwargs is None:
hdi_kwargs = {}
x = np.asarray(x)
x_shape = x.shape
if y is None and hdi_data is None:
raise ValueError("One of {y, hdi_data} is required")
if hdi_data is not None and y is not None:
warnings.warn("Both y and hdi_data arguments present, ignoring y")
elif hdi_data is not None:
hdi_prob = (
hdi_data.hdi.attrs.get("hdi_prob", np.nan) if hasattr(hdi_data, "hdi") else np.nan
)
if isinstance(hdi_data, Dataset):
data_vars = list(hdi_data.data_vars)
if len(data_vars) != 1:
raise ValueError(
"Found several variables in hdi_data. Only single variable Datasets are "
"supported."
)
hdi_data = hdi_data[data_vars[0]]
else:
y = np.asarray(y)
if hdi_prob is None:
hdi_prob = rcParams["stats.hdi_prob"]
elif not 1 >= hdi_prob > 0:
raise ValueError("The value of hdi_prob should be in the interval (0, 1]")
hdi_data = hdi(y, hdi_prob=hdi_prob, circular=circular, multimodal=False, **hdi_kwargs)
hdi_shape = hdi_data.shape
if hdi_shape[:-1] != x_shape:
msg = (
"Dimension mismatch for x: {} and hdi: {}. Check the dimensions of y and"
"hdi_kwargs to make sure they are compatible"
)
raise TypeError(msg.format(x_shape, hdi_shape))
if smooth:
if isinstance(x[0], np.datetime64):
raise TypeError("Cannot deal with x as type datetime. Recommend setting smooth=False.")
if smooth_kwargs is None:
smooth_kwargs = {}
smooth_kwargs.setdefault("window_length", 55)
smooth_kwargs.setdefault("polyorder", 2)
x_data = np.linspace(x.min(), x.max(), 200)
x_data[0] = (x_data[0] + x_data[1]) / 2
hdi_interp = griddata(x, hdi_data, x_data)
y_data = savgol_filter(hdi_interp, axis=0, **smooth_kwargs)
else:
idx = np.argsort(x)
x_data = x[idx]
y_data = hdi_data[idx]
hdiplot_kwargs = dict(
ax=ax,
x_data=x_data,
y_data=y_data,
color=color,
figsize=figsize,
plot_kwargs=plot_kwargs,
fill_kwargs=fill_kwargs,
backend_kwargs=backend_kwargs,
show=show,
)
if backend is None:
backend = rcParams["plot.backend"]
backend = backend.lower()
plot = get_plotting_function("plot_hdi", "hdiplot", backend)
ax = plot(**hdiplot_kwargs)
return ax
