From c66080768917cccc156faef3f2fd94d2021c8026 Mon Sep 17 00:00:00 2001
From: tomicapretto <tomicapretto@gmail.com>
Date: Tue, 7 Jul 2020 16:10:51 -0300
Subject: [PATCH 01/11] New KDE function. Added a new KDE function that can be
 used to estimate linear and circular data. Decided to put everything in a new
 file (instead of numeric_utils) because these functions required many
 utilities themselves. Also, updated all functions that used the old 
 function. The circular option is not yet passed to higher level functions.
 That is something to be done.

---
 arviz/kde_utils.py                            | 812 ++++++++++++++++++
 arviz/numeric_utils.py                        |  80 +-
 arviz/plots/backends/bokeh/bpvplot.py         |   9 +-
 arviz/plots/backends/bokeh/densityplot.py     |  11 +-
 arviz/plots/backends/bokeh/forestplot.py      |   7 +-
 arviz/plots/backends/bokeh/loopitplot.py      |   6 +-
 arviz/plots/backends/bokeh/ppcplot.py         |   6 +-
 arviz/plots/backends/bokeh/violinplot.py      |   4 +-
 arviz/plots/backends/matplotlib/bpvplot.py    |   9 +-
 .../plots/backends/matplotlib/densityplot.py  |  22 +-
 arviz/plots/backends/matplotlib/forestplot.py |   6 +-
 arviz/plots/backends/matplotlib/loopitplot.py |   4 +-
 arviz/plots/backends/matplotlib/ppcplot.py    |  18 +-
 arviz/plots/backends/matplotlib/violinplot.py |   6 +-
 arviz/plots/densityplot.py                    |   8 +-
 arviz/plots/distplot.py                       |   8 +-
 arviz/plots/energyplot.py                     |   9 +-
 arviz/plots/kdeplot.py                        |  23 +-
 arviz/plots/loopitplot.py                     |   5 +-
 arviz/plots/plot_utils.py                     |  18 +-
 arviz/plots/posteriorplot.py                  |  11 +-
 arviz/plots/violinplot.py                     |   7 +-
 arviz/stats/stats.py                          |   7 +-
 .../tests/base_tests/test_plots_matplotlib.py |  29 +-
 24 files changed, 927 insertions(+), 198 deletions(-)
 create mode 100644 arviz/kde_utils.py

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
new file mode 100644
index 0000000000..f83a6bc07f
--- /dev/null
+++ b/arviz/kde_utils.py
@@ -0,0 +1,812 @@
+"""Kernel density estimation functions for ArviZ."""
+
+import numpy as np
+from scipy.signal import gaussian, convolve
+from scipy.fftpack import fft
+from scipy.optimize import fsolve
+from scipy.optimize import brentq
+from scipy.special import ive
+import warnings
+
+# Own
+from .stats.stats_utils import histogram
+
+# Bandwidth functions ---------------------------------------------------------------
+# Linear KDE
+def bw_scott(x, x_std=None):
+    """
+    Scott's Rule
+    """
+    if x_std is None:
+        x_std = np.std(x)
+    bw = 1.06 * x_std * len(x) ** (-0.2)
+    return bw
+
+
+def bw_silverman(x, x_std=None):
+    """
+    Silverman's Rule.
+    """
+    if x_std is None:
+        x_std = np.std(x)
+    q75, q25 = np.percentile(x, [75, 25])
+    x_iqr = q75 - q25
+    a = min(x_std, x_iqr / 1.34)
+    bw = 0.9 * a * len(x) ** (-0.2)
+    return bw
+
+
+def bw_isj(x, grid_counts=None, x_range=None):
+    """
+    Improved Sheather and Jones method as explained in [1]
+    This is an internal version pretended to be used by the KDE estimator.
+    When used internally computation time is saved because
+    things like minimums, maximums and the grid are pre-computed.
+
+    References
+    ----------
+    .. [1] Kernel density estimation via diffusion.
+       Z. I. Botev, J. F. Grotowski, and D. P. Kroese.
+       Ann. Statist. 38 (2010), no. 5, 2916--2957.
+    """
+
+    x_len = len(x)
+    if x_range is None:
+        x_min = np.min(x)
+        x_max = np.max(x)
+        x_range = x_max - x_min
+
+    # Relative frequency per bin
+    if grid_counts is None:
+        x_std = np.std(x)
+        grid_len = 256
+        grid_min = x_min - 0.5 * x_std
+        grid_max = x_max + 0.5 * x_std
+        grid_counts, _, _ = histogram(x, grid_len, (grid_min, grid_max))
+    else:
+        grid_len = len(grid_counts) - 1
+
+    grid_relfreq = grid_counts / x_len
+
+    # Discrete cosine transform of the data
+    a_k = dct1d(grid_relfreq)
+
+    k_sq = np.arange(1, grid_len) ** 2
+    a_sq = a_k[range(1, grid_len)] ** 2
+
+    t = _root(_fixed_point, x_len, args=(x_len, k_sq, a_sq), x=x)
+    h = t ** 0.5 * x_range
+    return h
+
+def bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
+    """
+    Experimental bandwidth estimator.
+    """
+    return 0.5 * (bw_silverman(x, x_std) + bw_isj(x, grid_counts, x_range))
+
+
+def bw_circular(x):
+    """
+    Bandwidth selector for circular kernel density estimation
+    as introduced in [1].
+    This function implements a rule-of-thumb for choosing the bandwidth of
+    a von Mises kernel density estimator that assumes the underlying
+    distribution is von Mises.
+    It is analogous to Scott's rule for the Gaussian KDE.
+
+    Circular bandwidth has a different scale from linear bandwidth.
+    Unlike linear scale, low bandwidths are associated with oversmoothing
+    while high values are associated with undersmoothing.
+
+    References
+    ----------
+    [1] C.C Taylor (2008). Automatic bandwidth selection
+    for circular density estimation.
+    Computational Statistics and Data Analysis, 52, 7, 3493–3500.
+    """
+    x_len = len(x)
+    kappa = kappa_mle(x)
+    num = 3 * x_len * kappa ** 2 * ive(2, 2 * kappa)
+    den = 4 * np.pi ** 0.5 * ive(0, kappa) ** 2
+    return (num / den) ** 0.4
+
+
+BW_METHODS = {
+        "scott": bw_scott,
+        "silverman": bw_silverman,
+        "isj": bw_isj,
+        "experimental": bw_experimental
+    }
+
+
+def select_bw_method(method):
+    """
+    Selects a function to compute the bandwidth.
+    Also checks method `bw` is correctly specified.
+    Otherwise, throws an error.
+
+    Parameters
+    ----------
+    method : str
+        Method to estimate the bandwidth.
+
+    Returns
+    -------
+    bw_fun: function
+        Function to compute the bandwidth.
+    """
+
+    method_lower = method.lower()
+    if method_lower not in BW_METHODS.keys():
+        raise ValueError((
+            f"Unrecognized bandwidth method.\n"
+            f"Input is: {method}.\n"
+            f"Expected one of: {list(BW_METHODS.keys())}."
+        ))
+    bw_fun = BW_METHODS[method_lower]
+    return bw_fun
+
+
+def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
+    """
+    Computes bandwidth for a given data `x` and `bw`.
+    Also checks `bw` is correctly specified.
+
+    Parameters
+    ----------
+    x : 1-D numpy array
+        1 dimensional array of sample data from the
+        variable for which a density estimate is desired.
+    bw: int, float or str
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth.
+
+    Returns
+    -------
+    bw: float
+        Bandwidth
+    """
+    if isinstance(bw, bool):
+        raise ValueError((
+            f"`bw` must not be of type `bool`.\n"
+            f"Expected a positive numeric or one of the following strings:\n"
+            f"{list(BW_METHODS.keys())}."))
+    if isinstance(bw, (int, float)):
+        if bw < 0:
+            raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
+    elif isinstance(bw, str):
+        bw_fun = select_bw_method(bw)
+        bw_lower = bw.lower()
+        if bw_lower == "isj":
+            bw = bw_fun(x, grid_counts, x_range)
+        elif bw_lower in ["scott", "silverman"]:
+            bw = bw_fun(x, x_std)
+        elif bw_lower == "experimental":
+            bw = bw_fun(x, grid_counts, x_std, x_range)
+    else:
+        raise ValueError((
+            f"Unrecognized `bw` argument.\n"
+            f"Expected a positive numeric or one of the following strings:\n"
+            f"{list(BW_METHODS.keys())}."))
+    return bw
+
+
+# Misc utils
+def vonmises_pdf(x, mu, kappa):
+    assert kappa > 0, "Argument 'kappa' must be positive."
+    pdf = 1 / (2 * np.pi * ive(0, kappa)) * np.exp(np.cos(x - mu) - 1) ** kappa
+    return pdf
+
+
+def circular_mean(x):
+    sinr = np.sum(np.sin(x))
+    cosr = np.sum(np.cos(x))
+    mean = np.arctan2(sinr, cosr)
+    return mean
+
+
+def A1inv(x):
+    """
+    Inverse function of the ratio of the first and zeroth order
+    Bessel functions of the first kind.
+
+    Returns the value k, such that A1inv(x) = k, i.e. A1(k) = x.
+    """
+    if 0 <= x < 0.53:
+        return 2 * x + x ** 3 + (5 * x ** 5) / 6
+    elif x < 0.85:
+        return -0.4 + 1.39 * x + 0.43 / (1 - x)
+    else:
+        return 1 / (x ** 3 - 4 * x ** 2 + 3 * x)
+
+
+def kappa_mle(x):
+    mean = circular_mean(x)
+    kappa = A1inv(np.mean(np.cos(x - mean)))
+    return kappa
+
+
+def dct1d(x):
+    """
+    Discrete Cosine Transform in 1 Dimension
+
+    Parameters
+    ----------
+    x : numpy array
+        1 dimensional array of values for which the
+        DCT is desired
+
+    Returns
+    -------
+    output : DTC transformed values
+    """
+
+    x_len = len(x)
+
+    even_increasing = np.arange(0, x_len, 2)
+    odd_decreasing = np.arange(x_len - 1, 0, -2)
+
+    x = np.concatenate((x[even_increasing], x[odd_decreasing]))
+
+    w_1k = np.r_[
+        1,
+        (2 * np.exp(-(0 + 1j) * (np.arange(1, x_len)) * np.pi / (2 * x_len)))
+    ]
+    output = np.real(w_1k * fft(x))
+
+    return output
+
+
+def _fixed_point(t, N, k_sq, a_sq):
+    """
+    Implementation of the function t-zeta*gamma^[l](t) derived
+    from equation (30) in [1]
+
+    References
+    ----------
+    .. [1] Kernel density estimation via diffusion.
+       Z. I. Botev, J. F. Grotowski, and D. P. Kroese.
+       Ann. Statist. 38 (2010), no. 5, 2916--2957.
+    """
+
+    k_sq = np.asfarray(k_sq, dtype=np.float64)
+    a_sq = np.asfarray(a_sq, dtype=np.float64)
+
+    l = 7
+    f = np.sum(np.power(k_sq, l) * a_sq * np.exp(-k_sq * np.pi ** 2 * t))
+    f *= 0.5 * np.pi ** (2.0 * l)
+
+    for j in reversed(range(2, l)):
+        c1 = (1 + 0.5 ** (j + 0.5)) / 3
+        c2 = np.product(np.arange(1.0, 2 * j + 1, 2, dtype=np.float64))
+        c2 /= (np.pi / 2) ** 0.5
+        t_j = np.power((c1 * (c2 / (N * f))), (2.0 / (3.0 + 2.0 * j)))
+        f = np.sum(k_sq ** j * a_sq * np.exp(-k_sq * np.pi ** 2.0 * t_j))
+        f *= 0.5 * np.pi ** (2 * j)
+
+    out = t - (2 * N * np.pi ** 0.5 * f) ** (-0.4)
+    return out
+
+
+def _root(function, N, args, x):
+    # The idea here was borrowed from KDEpy implementation.
+    # The right bound is at most 0.01
+    found = 0
+    N = max(min(1050, N), 50)
+    tol = 10e-12 + 0.01 * (N - 50) / 1000
+
+    while found == 0:
+        try:
+            bw, res = brentq(function, 0, 0.01, args=args, full_output=True,
+                             disp=False)
+            found = 1 if res.converged else 0
+        except ValueError:
+            bw = 0
+            tol *= 2.0
+            found = 0
+        if bw <= 0:
+            warnings.warn(
+                "Improved Sheather-Jones did not converge to a positive value. "
+                "Using Silverman's rule instead.",
+                Warning
+            )
+            bw = (bw_silverman(x) / np.ptp(x)) ** 2
+            return bw
+        if tol >= 1:
+            warnings.warn(
+                "Improved Sheather-Jones did not converge. "
+                "Using Silverman's rule instead.",
+                Warning
+            )
+            bw = (bw_silverman(x) / np.ptp(x)) ** 2
+    return bw
+
+# KDE Utilities ---------------------------------------------------------------------
+
+def check_type(x):
+    """
+    Checks the input is of the correct type.
+    It only accepts numeric lists/numpy arrays of 1 dimension.
+    If input is not of the correct type, an informative message is thrown.
+
+    Parameters
+    ----------
+    x : Object whose type is checked before computing the KDE.
+
+    Returns
+    -------
+    x : 1-D numpy array
+        If no error is thrown, a 1 dimensional array of
+        sample data from the variable for which a density estimate is desired.
+
+    """
+    if not isinstance(x, (list, np.ndarray)):
+        raise ValueError((
+            f"`x` is of the wrong type.\n"
+            f"Can't produce a density estimator for {type(x)}.\n"
+            f"Please input a numeric list or numpy array."
+        ))
+
+    # Will raise an error if `x` can't be casted to numeric
+    x = np.asfarray(x)
+
+    x = x[np.isfinite(x)]
+    if x.size == 0:
+        raise ValueError("`x` does not contain any finite number.")
+
+    if x.ndim != 1:
+        raise ValueError((
+            f"Unsupported dimension number.\n"
+            f"Density estimator only works with 1-dimensional data, not {x.ndim}."
+        ))
+    return x
+
+
+def check_custom_lims(custom_lims, x_min, x_max):
+    """
+    Checks whether `custom_lims` are of the correct type.
+    It accepts numeric lists/tuples of length 2.
+
+    Parameters
+    ----------
+    custom_lims : Object whose type is checked.
+
+    Returns
+    -------
+    None: Object of type None
+
+    """
+    if not isinstance(custom_lims, (list, tuple)):
+        raise TypeError((
+            f"`custom_lims` must be a numeric list or tuple of length 2.\n"
+            f"Not an object of {type(custom_lims)}."
+        ))
+
+    if len(custom_lims) != 2:
+        raise AttributeError(
+            f"`len(custom_lims)` must be 2, not {len(custom_lims)}.")
+
+    any_bool = any(isinstance(i, bool) for i in custom_lims)
+    if any_bool:
+        raise TypeError(
+            "Elements of `custom_lims` must be numeric or None, not bool.")
+
+    if custom_lims[0] is None:
+        custom_lims[0] = x_min
+
+    if custom_lims[1] is None:
+        custom_lims[1] = x_max
+
+    types = (int, float, np.integer, np.float)
+    all_numeric = all(isinstance(i, types) for i in custom_lims)
+    if not all_numeric:
+        raise TypeError((
+            f"Elements of `custom_lims` must be numeric or None.\n"
+            f"At least one of them is not."
+        ))
+
+    if not custom_lims[0] < custom_lims[1]:
+        raise AttributeError(
+            f"`custom_lims[0]` must be smaller than `custom_lims[1]`.")
+
+    return custom_lims
+
+def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True,
+            bound_correction=False):
+    """
+    Computes the grid that bins the data used to estimate the density function
+
+    Parameters
+    ----------
+    x_min : float
+        Minimum value of the data
+    x_max: float
+        Maximum value of the data.
+    x_std: float
+        Standard deviation of the data.
+    extend_fct: bool
+        Indicates the factor by which `x_std` is multiplied
+        to extend the range of the data.
+    grid_len: int
+        Number of bins
+    custom_lims: tuple or list
+        Custom limits for the domain of the density estimation.
+        Must be numeric of length 2.
+    extend: bool, optional
+        Whether to extend the range of the data or not.
+        Default is True.
+    bound_correction: bool, optional
+        Whether the density estimations performs boundary correction or not.
+        This does not impacts directly in the output, but is used
+        to override `extend`.
+        Default is False.
+
+    Returns
+    -------
+    grid_len: int
+        Number of bins
+    grid_min: float
+        Minimum value of the grid
+    grid_max: float
+        Maximum value of the grid
+
+    """
+
+    # Set up number of bins. Should enable larger grids?
+    if grid_len > 4096:
+        grid_len = 4096
+    if grid_len < 100:
+        grid_len = 100
+    grid_len = int(grid_len)
+
+    # Set up domain
+    # `custom_lims` overrides `extend`
+    # `bound_correction` overrides `extend`
+    if custom_lims is not None:
+        custom_lims = check_custom_lims(custom_lims, x_min, x_max)
+        grid_min = custom_lims[0]
+        grid_max = custom_lims[1]
+    elif extend and not bound_correction:
+        grid_extend = extend_fct * x_std
+        grid_min = x_min - grid_extend
+        grid_max = x_max + grid_extend
+    else:
+        grid_min = x_min
+        grid_max = x_max
+    return grid_min, grid_max, grid_len
+
+# KDE Functions --------------------------------------------------------------------
+
+def kde(x, circular=False, **kwargs):
+    """
+    1 dimensional density estimation.
+    It is a wrapper around `kde_linear()` and `kde_circular()`.
+
+    Parameters
+    ----------
+    x : 1D numpy array
+        Data used to calculate the density estimation.
+        Theoritically it is a random sample obtained from $f$,
+        the true probability density function we aim to estimate.
+    circular: bool, optional
+        Whether `x` is a circular variable or not. Defaults to False.
+    **kwargs: Arguments passed to `kde_linear()` and `kde_circular()`.
+        See their documentation for more info.
+
+    Returns
+    -------
+    grid : Gridded numpy array for the x values.
+    pdf : Numpy array for the density estimates.
+    bw: optional, the estimated bandwidth.
+    """
+    if circular:
+        kde_fun = kde_circular
+    else:
+        kde_fun = kde_linear
+
+    return kde_fun(x, **kwargs)
+
+
+def kde_linear(
+    x,
+    bw="experimental",
+    adaptive=False,
+    extend=False,
+    bound_correction=True,
+    extend_fct=0,
+    bw_fct=1,
+    bw_return=False,
+    custom_lims=None,
+    cumulative=False,
+    grid_len=512
+):
+    """
+    1 dimensional density estimation for linear data.
+
+    Given an array of data points `x` it returns an estimate of
+    the probability density function that generated the samples in `x`.
+
+    Parameters
+    ----------
+    x : 1D numpy array
+        Data used to calculate the density estimation.
+        Theoritically it is a random sample obtained from $f$,
+        the true probability density function we aim to estimate.
+    bw: int, float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental".
+        Defaults to "experimental".
+    adaptive: boolean, optional
+        Indicates if the bandwidth is adaptative or not.
+        It is the recommended approach when there are multiple modalities
+        with different spread.
+        It is not compatible with convolution. Defaults to False.
+    extend: boolean, optional
+        Whether to extend the observed range for `x` in the estimation.
+        It extends each bound by a multiple of the standard deviation of `x`
+        given by `extend_fct`. Defaults to False.
+    bound_correction: boolean, optional
+        Whether to perform boundary correction on the bounds of `x` or not.
+        Defaults to True.
+    extend_fct: float, optional
+        Number of standard deviations used to widen the
+        lower and upper bounds of `x`. Defaults to 0.5.
+    bw_fct: float, optional
+        A value that multiplies `bw` which enables tuning smoothness by hand.
+        Must be positive. Values below 1 decrease smoothness while values
+        above 1 decrease it. Defaults to 1 (no modification).
+    bw_return: bool, optional
+        Whether to return the estimated bandwidth in addition to the
+        other objects. Defaults to False.
+    custom_lims: list or tuple, optional
+        A list or tuple of length 2 indicating custom bounds
+        for the range of `x`. Defaults to None which disables custom bounds.
+    cumulative: bool, optional
+        Whether return the PDF or the cumulative PDF. Defaults to False.
+    grid_len: int, optional
+        The number of intervals used to bin the data points
+        (a.k.a. the length of the grid used in the estimation)
+        Defaults to 512.
+
+    Returns
+    -------
+    grid : Gridded numpy array for the x values.
+    pdf : Numpy array for the density estimates.
+    bw: optional, the estimated bandwidth.
+    """
+
+    # Check `x` is from appropiate type
+    x = check_type(x)
+
+    # Assert `bw_fct` is numeric and positive
+    assert isinstance(bw_fct, (int, float))
+    assert bw_fct > 0
+
+    # Preliminary calculations
+    x_len = len(x)
+    x_min = x.min()
+    x_max = x.max()
+    x_std = (((x ** 2).sum() / x_len) - (x.sum() / x_len) ** 2) ** 0.5
+    x_range = x_max - x_min
+
+    # Determine grid
+    grid_min, grid_max, grid_len = get_grid(
+        x_min, x_max, x_std, extend_fct, grid_len,
+        custom_lims, extend, bound_correction
+    )
+    grid_counts, _, grid_edges = histogram(x, grid_len, (grid_min, grid_max))
+
+    # Bandwidth estimation
+    bw = bw_fct * get_bw(x, bw, grid_counts, x_std, x_range)
+
+    # Density estimation
+    if adaptive:
+        grid, pdf = _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len,
+                                  bound_correction)
+    else:
+        grid, pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len,
+                                     bound_correction)
+
+    if cumulative:
+        pdf = pdf.cumsum() / pdf.sum()
+
+    if bw_return:
+        return grid, pdf, bw
+    else:
+        return grid, pdf
+
+
+def kde_circular(
+    x,
+    bw="default",
+    bw_fct=1,
+    bw_return=False,
+    custom_lims=None,
+    cumulative=False,
+    grid_len=512
+):
+    """
+    1 dimensional density estimation for circular data.
+
+    Given an array of data points `x` measured in radians,
+    it returns an estimate of the probability density function that generated
+    the samples in `x`.
+
+    Parameters
+    ----------
+    x : 1D numpy array
+        Data used to calculate the density estimation.
+        Theoritically it is a random sample obtained from $f$,
+        the true probability density function we aim to estimate.
+    bw: int, float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one "default". Defaults to "default".
+    bw_fct: float, optional
+        A value that multiplies `bw` which enables tuning smoothness by hand.
+        Must be positive. Values above 1 decrease smoothness while values
+        below 1 decrease it. Defaults to 1 (no modification).
+    bw_return: bool, optional
+        Whether to return the estimated bandwidth in addition to the
+        other objects. Defaults to False.
+    custom_lims: list or tuple, optional
+        A list or tuple of length 2 indicating custom bounds
+        for the range of `x`. Defaults to None which means the estimation
+        limits are [-pi, pi].
+    cumulative: bool, optional
+        Whether return the PDF or the cumulative PDF. Defaults to False.
+    grid_len: int, optional
+        The number of intervals used to bin the data points
+        (a.k.a. the length of the grid used in the estimation)
+        Defaults to 512.
+    """
+
+    x = check_type(x)
+
+    # All values between -pi and pi
+    x[x > np.pi] = x[x > np.pi] - 2 * np.pi
+    x[x < -np.pi] = x[x < -np.pi] + 2 * np.pi
+
+    # Assert `bw_fct` is numeric and positive
+    assert isinstance(bw_fct, (int, float))
+    assert bw_fct > 0
+
+    # Determine bandwidth
+    if isinstance(bw, bool):
+        raise ValueError((
+            "`bw` can't be of type `bool`.\n"
+            "Expected a positive numeric or 'default'"
+        ))
+    if isinstance(bw, (int, float)):
+        if bw < 0:
+            raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
+    if isinstance(bw, str):
+        if bw == "default":
+            bw = bw_circular(x)
+        else:
+            raise ValueError((
+                f"`bw` must be a positive numeric or `default`, not {bw}"
+            ))
+    bw *= bw_fct
+
+    # Determine grid
+    if custom_lims is not None:
+        custom_lims = check_custom_lims(custom_lims, x.min(), x.max())
+        grid_min = custom_lims[0]
+        grid_max = custom_lims[1]
+        assert grid_min >= -np.pi, "Lower limit can't be smaller than -pi"
+        assert grid_max <= np.pi, "Upper limit can't be larger than pi"
+    else:
+        grid_min = -np.pi
+        grid_max = np.pi
+
+    bins = np.linspace(grid_min, grid_max, grid_len + 1)
+    bin_counts, _, bin_edges = histogram(x, bins=bins)
+    grid = 0.5 * (bin_edges[1:] + bin_edges[:-1])
+
+    kern = vonmises_pdf(x=grid, mu=0, kappa=bw)
+    pdf = np.fft.fftshift(np.fft.irfft(np.fft.rfft(kern) * np.fft.rfft(bin_counts)))
+    pdf /= len(x)
+
+    if cumulative:
+        pdf = pdf.cumsum() / pdf.sum()
+
+    if bw_return:
+        return grid, pdf, bw
+    else:
+        return grid, pdf
+
+
+def _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
+    """
+    1 dimensional Gaussian kernel density estimation via
+    convolution of the binned relative frequencies and a Gaussian filter.
+    This is an internal function used by `kde()`.
+    """
+
+    # Calculate relative frequencies per bin
+    bin_width = grid_edges[1] - grid_edges[0]
+    f = grid_counts / bin_width / len(x)
+
+    # Bandwidth must consider the bin width
+    bw /= bin_width
+
+    # See: https://stackoverflow.com/questions/2773606/gaussian-filter-in-matlab
+    kernel_n = int(bw * 2 * np.pi)
+
+    # Temporal fix?
+    if kernel_n == 0:
+        kernel_n = 1
+    kernel = gaussian(kernel_n, bw)
+
+    if bound_correction:
+        npad = int(grid_len / 5)
+        f = np.concatenate([f[npad - 1::-1], f, f[grid_len:grid_len - npad - 1:-1]])
+        pdf = convolve(f, kernel, mode="same", method="direct")[npad:npad + grid_len]
+        pdf /= bw * (2 * np.pi) ** 0.5
+    else:
+        pdf = convolve(f, kernel, mode="same", method="direct")
+        pdf /= (bw * (2 * np.pi) ** 0.5)
+
+    grid = (grid_edges[1:] + grid_edges[:-1]) / 2
+    return grid, pdf
+
+
+def _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
+    """
+    1 dimensional adaptive Gaussian kernel density estimation.
+    The implementation uses the binning technique.
+    Since there is not an unique `bw`, the convolution is not possible.
+    The alternative implemented in this function is known as Abramson's method.
+    This is an internal function used by `kde()`.
+    """
+    # Pilot computations used for bandwidth adjustment
+    pilot_grid, pilot_pdf = _kde_convolution(x, bw, grid_edges, grid_counts,
+                                             grid_len, bound_correction)
+
+    # Adds to avoid np.log(0) and zero division
+    pilot_pdf += 1e-9
+
+    # Determine the modification factors
+    pdf_interp = np.interp(x, pilot_grid, pilot_pdf)
+    geom_mean = np.exp(np.mean(np.log(pdf_interp)))
+
+    # Power of c = 0.5 -> Abramson's method
+    adj_factor = (geom_mean / pilot_pdf) ** 0.5
+    bw_adj = bw * adj_factor
+
+    # Estimation of Gaussian KDE via binned method (convolution not possible)
+    grid = pilot_grid
+
+    if bound_correction:
+        grid_npad = int(grid_len / 5)
+        grid_width = grid_edges[1] - grid_edges[0]
+        grid_pad = grid_npad * grid_width
+        grid_padded = np.linspace(
+            grid_edges[0] - grid_pad,
+            grid_edges[grid_len - 1] + grid_pad,
+            num = grid_len + 2 * grid_npad
+        )
+        grid_counts = np.concatenate([
+            grid_counts[grid_npad - 1:: -1],
+            grid_counts,
+            grid_counts[grid_len:grid_len - grid_npad - 1: -1]]
+        )
+        bw_adj = np.concatenate([
+            bw_adj[grid_npad - 1:: -1],
+            bw_adj,
+            bw_adj[grid_len:grid_len - grid_npad - 1: -1]]
+        )
+        pdf_mat = ((grid_padded - grid_padded[:, None]) / bw_adj[:, None])
+        pdf_mat = np.exp(-0.5 * pdf_mat ** 2) * grid_counts[:, None]
+        pdf_mat /= ((2 * np.pi) ** 0.5 * bw_adj[:, None])
+        pdf = np.sum(pdf_mat[:, grid_npad:grid_npad + grid_len], axis=0) / len(x)
+
+    else:
+        pdf_mat = ((grid - grid[:, None]) / bw_adj[:, None])
+        pdf_mat = np.exp(-0.5 * pdf_mat ** 2) * grid_counts[:, None]
+        pdf_mat /= ((2 * np.pi) ** 0.5 * bw_adj[:, None])
+        pdf = np.sum(pdf_mat, axis=0) / len(x)
+
+    return grid, pdf
diff --git a/arviz/numeric_utils.py b/arviz/numeric_utils.py
index c849b251d0..35f66266e5 100644
--- a/arviz/numeric_utils.py
+++ b/arviz/numeric_utils.py
@@ -2,84 +2,11 @@
 import warnings
 import numpy as np
 from scipy.signal import convolve, convolve2d
-from scipy.signal.windows import gaussian
 from scipy.sparse import coo_matrix
 
 from .stats.stats_utils import histogram
 from .utils import _stack, _dot, _cov
 
-
-def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):
-    """Fast Fourier transform-based Gaussian kernel density estimate (KDE).
-
-    The code was adapted from https://github.com/mfouesneau/faststats
-
-    Parameters
-    ----------
-    x : Numpy array or list
-    cumulative : bool
-        If true, estimate the cdf instead of the pdf
-    bw : float
-        Bandwidth scaling factor for the KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule
-        of thumb (the default rule used by SciPy).
-    xmin : float
-        Manually set lower limit.
-    xmax : float
-        Manually set upper limit.
-
-    Returns
-    -------
-    density: A gridded 1D KDE of the input points (x)
-    xmin: minimum value of x
-    xmax: maximum value of x
-    """
-    x = np.asarray(x, dtype=float)
-    x = x[np.isfinite(x)]
-    if x.size == 0:
-        warnings.warn("kde plot failed, you may want to check your data")
-        return np.array([np.nan]), np.nan, np.nan
-
-    len_x = len(x)
-    n_points = 200 if (xmin or xmax) is None else 500
-
-    if xmin is None:
-        xmin = np.min(x)
-    if xmax is None:
-        xmax = np.max(x)
-
-    assert np.min(x) >= xmin
-    assert np.max(x) <= xmax
-
-    log_len_x = np.log(len_x) * bw
-
-    n_bins = min(int(len_x ** (1 / 3) * log_len_x * 2), n_points)
-    if n_bins < 2:
-        warnings.warn("kde plot failed, you may want to check your data")
-        return np.array([np.nan]), np.nan, np.nan
-
-    # hist, bin_edges = np.histogram(x, bins=n_bins, range=(xmin, xmax))
-    # grid = hist / (hist.sum() * np.diff(bin_edges))
-
-    _, grid, _ = histogram(x, n_bins, range_hist=(xmin, xmax))
-
-    scotts_factor = len_x ** (-0.2)
-    kern_nx = int(scotts_factor * 2 * np.pi * log_len_x)
-    kernel = gaussian(kern_nx, scotts_factor * log_len_x)
-
-    npad = min(n_bins, 2 * kern_nx)
-    grid = np.concatenate([grid[npad:0:-1], grid, grid[n_bins : n_bins - npad : -1]])
-    density = convolve(grid, kernel, mode="same", method="direct")[npad : npad + n_bins]
-    norm_factor = (2 * np.pi * log_len_x ** 2 * scotts_factor ** 2) ** 0.5
-
-    density /= norm_factor
-
-    if cumulative:
-        density = density.cumsum() / density.sum()
-
-    return density, xmin, xmax
-
-
 def _fast_kde_2d(x, y, gridsize=(128, 128), circular=False):
     """
     2D fft-based Gaussian kernel density estimate (KDE).
@@ -185,7 +112,8 @@ def get_bins(values):
     bins_sturges = (x_max - x_min) / (np.log2(values.size) + 1)
 
     # The Freedman-Diaconis histogram bin estimator.
-    iqr = np.subtract(*np.percentile(values, [75, 25]))  # pylint: disable=assignment-from-no-return
+    iqr = np.subtract(
+        *np.percentile(values, [75, 25]))  # pylint: disable=assignment-from-no-return
     bins_fd = 2 * iqr * values.size ** (-1 / 3)
 
     width = np.round(np.max([1, bins_sturges, bins_fd])).astype(int)
@@ -202,12 +130,12 @@ def _sturges_formula(dataset, mult=1):
     Parameters
     ----------
     dataset: xarray.DataSet
-        Must have the `draw` dimension
+       Must have the `draw` dimension
 
     mult: float
         Used to scale the number of bins up or down. Default is 1 for Sturges' formula.
 
-    Returns
+     Returns
     -------
     int
         Number of bins to use
diff --git a/arviz/plots/backends/bokeh/bpvplot.py b/arviz/plots/backends/bokeh/bpvplot.py
index bcb5c8ac28..d40560fd95 100644
--- a/arviz/plots/backends/bokeh/bpvplot.py
+++ b/arviz/plots/backends/bokeh/bpvplot.py
@@ -14,8 +14,7 @@
     is_valid_quantile,
     vectorized_to_hex,
 )
-from ....numeric_utils import _fast_kde
-
+from ....kde_utils import kde
 
 def plot_bpv(
     ax,
@@ -91,8 +90,7 @@ def plot_bpv(
 
         if kind == "p_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=-1)
-            tstat_pit_dens, xmin, xmax = _fast_kde(tstat_pit)
-            x_s = np.linspace(xmin, xmax, len(tstat_pit_dens))
+            x_s, tstat_pit_dens = kde(tstat_pit)
             ax_i.line(x_s, tstat_pit_dens, line_width=linewidth, line_color=color)
             # ax_i.set_yticks([])
             if reference is not None:
@@ -112,8 +110,7 @@ def plot_bpv(
 
         elif kind == "u_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=0)
-            tstat_pit_dens, xmin, xmax = _fast_kde(tstat_pit)
-            x_s = np.linspace(xmin, xmax, len(tstat_pit_dens))
+            x_s, tstat_pit_dens = kde(tstat_pit)
             ax_i.line(x_s, tstat_pit_dens, line_color=color)
             if reference is not None:
                 if reference == "analytical":
diff --git a/arviz/plots/backends/bokeh/densityplot.py b/arviz/plots/backends/bokeh/densityplot.py
index a5c53a7757..6a50881e5a 100644
--- a/arviz/plots/backends/bokeh/densityplot.py
+++ b/arviz/plots/backends/bokeh/densityplot.py
@@ -16,8 +16,8 @@
     vectorized_to_hex,
 )
 from ....stats import hdi
-from ....numeric_utils import _fast_kde, histogram, get_bins
-
+from ....numeric_utils import histogram, get_bins
+from ....kde_utils import kde
 
 def plot_density(
     ax,
@@ -144,11 +144,10 @@ def _d_helper(
         else:
             new_vec = vec
 
-        density, xmin, xmax = _fast_kde(new_vec, bw=bw)
+        x, density = kde(new_vec, bw_fct=bw)
         density *= hdi_prob
-        x = np.linspace(xmin, xmax, len(density))
-        ymin = density[0]
-        ymax = density[-1]
+        xmin, xmax = x[0], x[-1]
+        ymin, ymax = density[0], density[-1]
 
         if outline:
             plotted.append(ax.line(x, density, line_color=color, line_width=line_width, **extra))
diff --git a/arviz/plots/backends/bokeh/forestplot.py b/arviz/plots/backends/bokeh/forestplot.py
index 79b5cec338..9e151011d5 100644
--- a/arviz/plots/backends/bokeh/forestplot.py
+++ b/arviz/plots/backends/bokeh/forestplot.py
@@ -16,10 +16,12 @@
 from ....rcparams import rcParams
 from ....stats import hdi
 from ....stats.diagnostics import _ess, _rhat
-from ....numeric_utils import _fast_kde, histogram, get_bins
+from ....numeric_utils import histogram, get_bins
+from ....kde_utils import kde
 from ....utils import conditional_jit
 
 
+
 def pairwise(iterable):
     """From itertools cookbook. [a, b, c, ...] -> (a, b), (b, c), ..."""
     first, second = tee(iterable)
@@ -570,9 +572,8 @@ def ridgeplot(self, mult, ridgeplot_kind):
                 _, density, x = histogram(values, bins=bins)
                 x = x[:-1]
             elif kind == "density":
-                density, lower, upper = _fast_kde(values)
+                x, density = kde(values)
                 density_q = density.cumsum() / density.sum()
-                x = np.linspace(lower, upper, len(density))
 
             xvals.append(x)
             pdfs.append(density)
diff --git a/arviz/plots/backends/bokeh/loopitplot.py b/arviz/plots/backends/bokeh/loopitplot.py
index d198822d95..5eeec57d75 100644
--- a/arviz/plots/backends/bokeh/loopitplot.py
+++ b/arviz/plots/backends/bokeh/loopitplot.py
@@ -7,10 +7,9 @@
 
 from . import backend_kwarg_defaults
 from .. import show_layout
-from ...kdeplot import _fast_kde
+from ....kde_utils import kde
 from ...plot_utils import _scale_fig_size
 
-
 def plot_loo_pit(
     ax,
     figsize,
@@ -184,8 +183,7 @@ def plot_loo_pit(
             )
         else:
             for idx in range(n_unif):
-                unif_density, xmin, xmax = _fast_kde(unif[idx, :])
-                x_s = np.linspace(xmin, xmax, len(unif_density))
+                x_s, unif_density = kde(unif[idx, :])
                 ax.line(
                     x_s,
                     unif_density,
diff --git a/arviz/plots/backends/bokeh/ppcplot.py b/arviz/plots/backends/bokeh/ppcplot.py
index c4e9461e6e..34b58b31ed 100644
--- a/arviz/plots/backends/bokeh/ppcplot.py
+++ b/arviz/plots/backends/bokeh/ppcplot.py
@@ -3,9 +3,10 @@
 
 from . import backend_kwarg_defaults
 from .. import show_layout
-from ...kdeplot import plot_kde, _fast_kde
+from ...kdeplot import plot_kde
 from ...plot_utils import _create_axes_grid, _scale_fig_size
 from ....numeric_utils import histogram, get_bins
+from ....kde_utils import kde
 
 
 def plot_ppc(
@@ -93,8 +94,7 @@ def plot_ppc(
             for vals in pp_sampled_vals:
                 vals = np.array([vals]).flatten()
                 if dtype == "f":
-                    pp_density, lower, upper = _fast_kde(vals)
-                    pp_x = np.linspace(lower, upper, len(pp_density))
+                    pp_x, pp_density = kde(vals)
                     pp_densities.append(pp_density)
                     pp_xs.append(pp_x)
                 else:
diff --git a/arviz/plots/backends/bokeh/violinplot.py b/arviz/plots/backends/bokeh/violinplot.py
index e13f40cc3a..fc58c1083f 100644
--- a/arviz/plots/backends/bokeh/violinplot.py
+++ b/arviz/plots/backends/bokeh/violinplot.py
@@ -4,6 +4,7 @@
 
 from . import backend_kwarg_defaults
 from .. import show_layout
+from ....kde_utils import kde
 from ....numeric_utils import _fast_kde, histogram, get_bins
 from ...plot_utils import make_label, _create_axes_grid, _scale_fig_size
 from ....stats import hdi
@@ -101,8 +102,7 @@ def plot_violin(
 
 def _violinplot(val, rug, shade, bw, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    density, low_b, up_b = _fast_kde(val, bw=bw)
-    x = np.linspace(low_b, up_b, len(density))
+    x, density = kde(val, bw_fct=bw)
 
     if not rug:
         x = np.concatenate([x, x[::-1]])
diff --git a/arviz/plots/backends/matplotlib/bpvplot.py b/arviz/plots/backends/matplotlib/bpvplot.py
index d02824a06b..e2b4c44e94 100644
--- a/arviz/plots/backends/matplotlib/bpvplot.py
+++ b/arviz/plots/backends/matplotlib/bpvplot.py
@@ -13,8 +13,7 @@
     is_valid_quantile,
     matplotlib_kwarg_dealiaser,
 )
-from ....numeric_utils import _fast_kde
-
+from ....kde_utils import kde
 
 def plot_bpv(
     ax,
@@ -77,8 +76,7 @@ def plot_bpv(
 
         if kind == "p_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=-1)
-            tstat_pit_dens, xmin, xmax = _fast_kde(tstat_pit)
-            x_s = np.linspace(xmin, xmax, len(tstat_pit_dens))
+            x_s, tstat_pit_dens = kde(tstat_pit)
             ax_i.plot(x_s, tstat_pit_dens, linewidth=linewidth, color=color)
             ax_i.set_yticks([])
             if reference is not None:
@@ -97,8 +95,7 @@ def plot_bpv(
 
         elif kind == "u_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=0)
-            tstat_pit_dens, xmin, xmax = _fast_kde(tstat_pit)
-            x_s = np.linspace(xmin, xmax, len(tstat_pit_dens))
+            x_s, tstat_pit_dens = kde(tstat_pit)
             ax_i.plot(x_s, tstat_pit_dens, color=color)
             if reference is not None:
                 if reference == "analytical":
diff --git a/arviz/plots/backends/matplotlib/densityplot.py b/arviz/plots/backends/matplotlib/densityplot.py
index 843078e044..e303046978 100644
--- a/arviz/plots/backends/matplotlib/densityplot.py
+++ b/arviz/plots/backends/matplotlib/densityplot.py
@@ -12,7 +12,8 @@
     _scale_fig_size,
     calculate_point_estimate,
 )
-from ....numeric_utils import _fast_kde, get_bins
+from ....numeric_utils import get_bins
+from ....kde_utils import kde
 
 
 def plot_density(
@@ -101,7 +102,7 @@ def _d_helper(
     vec,
     vname,
     color,
-    bw,
+    bw_fct,
     titlesize,
     xt_labelsize,
     linewidth,
@@ -123,10 +124,10 @@ def _d_helper(
         variable name
     color : str
         matplotlib color
-    bw : float
-        Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the
-        KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb
-        (the default used rule by SciPy).
+    bw_fct: float, optional
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
     titlesize : float
         font size for title
     xt_labelsize : float
@@ -152,11 +153,10 @@ def _d_helper(
         else:
             new_vec = vec
 
-        density, xmin, xmax = _fast_kde(new_vec, bw=bw)
+        x, density = kde(new_vec, bw_fct=bw_fct)
         density *= hdi_prob
-        x = np.linspace(xmin, xmax, len(density))
-        ymin = density[0]
-        ymax = density[-1]
+        xmin, xmax = x[0], x[-1]
+        ymin, ymax = density[0], density[-1]
 
         if outline:
             ax.plot(x, density, color=color, lw=linewidth)
@@ -179,7 +179,7 @@ def _d_helper(
         ax.plot(xmax, 0, hdi_markers, color=color, markeredgecolor="k", markersize=markersize)
 
     if point_estimate is not None:
-        est = calculate_point_estimate(point_estimate, vec, bw)
+        est = calculate_point_estimate(point_estimate, vec, bw_fct)
         ax.plot(est, 0, "o", color=color, markeredgecolor="k", markersize=markersize)
 
     ax.set_yticks([])
diff --git a/arviz/plots/backends/matplotlib/forestplot.py b/arviz/plots/backends/matplotlib/forestplot.py
index 2f0aabf445..a2afaed95b 100644
--- a/arviz/plots/backends/matplotlib/forestplot.py
+++ b/arviz/plots/backends/matplotlib/forestplot.py
@@ -9,7 +9,8 @@
 from . import backend_kwarg_defaults, backend_show
 from ....stats import hdi
 from ....stats.diagnostics import _ess, _rhat
-from ....numeric_utils import _fast_kde, histogram, get_bins
+from ....numeric_utils import histogram, get_bins
+from ....kde_utils import kde
 from ...plot_utils import _scale_fig_size, xarray_var_iter, make_label
 from ....utils import conditional_jit
 
@@ -518,9 +519,8 @@ def ridgeplot(self, mult, ridgeplot_kind):
                 density_q = density.cumsum() / density.sum()
                 x = x[:-1]
             elif kind == "density":
-                density, lower, upper = _fast_kde(values)
+                x, density = kde(values)
                 density_q = density.cumsum() / density.sum()
-                x = np.linspace(lower, upper, len(density))
 
             xvals.append(x)
             pdfs.append(density)
diff --git a/arviz/plots/backends/matplotlib/loopitplot.py b/arviz/plots/backends/matplotlib/loopitplot.py
index 8ef0568238..9e62ec34f2 100644
--- a/arviz/plots/backends/matplotlib/loopitplot.py
+++ b/arviz/plots/backends/matplotlib/loopitplot.py
@@ -10,6 +10,7 @@
     _scale_fig_size,
     matplotlib_kwarg_dealiaser,
 )
+from ....kde_utils import kde
 from ....numeric_utils import _fast_kde
 
 
@@ -118,8 +119,7 @@ def plot_loo_pit(
             ax.axhspan(*hdi_odds, **hdi_kwargs)
         else:
             for idx in range(n_unif):
-                unif_density, xmin, xmax = _fast_kde(unif[idx, :])
-                x_s = np.linspace(xmin, xmax, len(unif_density))
+                x_s, unif_density = kde(unif[idx, :])
                 x_ss[idx] = x_s
                 u_dens[idx] = unif_density
             ax.plot(x_ss.T, u_dens.T, **plot_unif_kwargs)
diff --git a/arviz/plots/backends/matplotlib/ppcplot.py b/arviz/plots/backends/matplotlib/ppcplot.py
index 8ebdfc40ce..abedc3d7d3 100644
--- a/arviz/plots/backends/matplotlib/ppcplot.py
+++ b/arviz/plots/backends/matplotlib/ppcplot.py
@@ -12,8 +12,8 @@
     _create_axes_grid,
     _scale_fig_size,
 )
-from ....numeric_utils import _fast_kde, histogram, get_bins
-
+from ....numeric_utils import histogram, get_bins
+from ....kde_utils import kde
 _log = logging.getLogger(__name__)
 
 
@@ -149,8 +149,7 @@ def plot_ppc(
             for vals in pp_sampled_vals:
                 vals = np.array([vals]).flatten()
                 if dtype == "f":
-                    pp_density, lower, upper = _fast_kde(vals)
-                    pp_x = np.linspace(lower, upper, len(pp_density))
+                    pp_x, pp_density = kde(vals)
                     pp_densities.append(pp_density)
                     pp_xs.append(pp_x)
                 else:
@@ -368,18 +367,13 @@ def _set_animation(
     if kind == "kde":
         length = len(pp_sampled_vals)
         if dtype == "f":
-            y_vals, lower, upper = _fast_kde(pp_sampled_vals[0])
-            x_vals = np.linspace(lower, upper, len(y_vals))
-
-            max_max = max([max(_fast_kde(pp_sampled_vals[i])[0]) for i in range(length)])
-
+            x_vals, y_vals = kde(pp_sampled_vals[0])
+            max_max = max([max(kde(pp_sampled_vals[i])[1]) for i in range(length)])
             ax.set_ylim(0, max_max)
-
             (line,) = ax.plot(x_vals, y_vals, **plot_kwargs)
 
             def animate(i):
-                y_vals, lower, upper = _fast_kde(pp_sampled_vals[i])
-                x_vals = np.linspace(lower, upper, len(y_vals))
+                x_vals, y_vals = kde(pp_sampled_vals[i])
                 line.set_data(x_vals, y_vals)
                 return line
 
diff --git a/arviz/plots/backends/matplotlib/violinplot.py b/arviz/plots/backends/matplotlib/violinplot.py
index 78ac2620b7..758aad9af4 100644
--- a/arviz/plots/backends/matplotlib/violinplot.py
+++ b/arviz/plots/backends/matplotlib/violinplot.py
@@ -4,7 +4,8 @@
 
 from . import backend_show
 from ....stats import hdi
-from ....numeric_utils import _fast_kde, histogram, get_bins
+from ....kde_utils import kde
+from ....numeric_utils import histogram, get_bins
 from ...plot_utils import make_label, _create_axes_grid, matplotlib_kwarg_dealiaser, _scale_fig_size
 
 
@@ -86,8 +87,7 @@ def plot_violin(
 
 def _violinplot(val, rug, shade, bw, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    density, low_b, up_b = _fast_kde(val, bw=bw)
-    x = np.linspace(low_b, up_b, len(density))
+    x, density = kde(val, bw_fct=bw)
 
     if not rug:
         x = np.concatenate([x, x[::-1]])
diff --git a/arviz/plots/densityplot.py b/arviz/plots/densityplot.py
index 8c4e1a1fb4..2955e08745 100644
--- a/arviz/plots/densityplot.py
+++ b/arviz/plots/densityplot.py
@@ -25,7 +25,7 @@ def plot_density(
     outline=True,
     hdi_markers="",
     shade=0.0,
-    bw=4.5,
+    bw=1,
     figsize=None,
     textsize=None,
     ax=None,
@@ -78,9 +78,9 @@ def plot_density(
         Alpha blending value for the shaded area under the curve, between 0 (no shade) and 1
         (opaque). Defaults to 0.
     bw : Optional[float]
-        Bandwidth scaling factor for the KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule
-        of thumb (the default rule used by SciPy).
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
     figsize : Optional[Tuple[int, int]]
         Figure size. If None it will be defined automatically.
     textsize: Optional[float]
diff --git a/arviz/plots/distplot.py b/arviz/plots/distplot.py
index bf222810d2..707a4ffb2b 100644
--- a/arviz/plots/distplot.py
+++ b/arviz/plots/distplot.py
@@ -15,7 +15,7 @@ def plot_dist(
     label=None,
     rotated=False,
     rug=False,
-    bw=4.5,
+    bw=1,
     quantiles=None,
     contour=True,
     fill_last=True,
@@ -59,9 +59,9 @@ def plot_dist(
     rug : bool
         If True adds a rugplot. Defaults to False. Ignored for 2D KDE
     bw : float
-        Bandwidth scaling factor for 1D KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's
-        rule of thumb (the default rule used by SciPy).
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
     quantiles : list
         Quantiles in ascending order used to segment the KDE. Use [.25, .5, .75] for quartiles.
         Defaults to None.
diff --git a/arviz/plots/energyplot.py b/arviz/plots/energyplot.py
index c607d19086..16df90e836 100644
--- a/arviz/plots/energyplot.py
+++ b/arviz/plots/energyplot.py
@@ -12,7 +12,7 @@ def plot_energy(
     legend=True,
     fill_alpha=(1, 0.75),
     fill_color=("C0", "C5"),
-    bw=4.5,
+    bw=1,
     textsize=None,
     fill_kwargs=None,
     plot_kwargs=None,
@@ -44,9 +44,10 @@ def plot_energy(
         Color for Marginal energy distribution and Energy transition distribution.
         Defaults to ('C0', 'C5')
     bw : float
-        Bandwidth scaling factor for the KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule
-        of thumb (the default rule used by SciPy). Only works if `kind='kde'`
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
+        Only works if `kind='kde'`
     textsize: float
         Text size scaling factor for labels, titles and lines. If None it will be autoscaled based
         on figsize.
diff --git a/arviz/plots/kdeplot.py b/arviz/plots/kdeplot.py
index c883c36247..7bfd21870b 100644
--- a/arviz/plots/kdeplot.py
+++ b/arviz/plots/kdeplot.py
@@ -3,7 +3,8 @@
 import xarray as xr
 
 from ..data import InferenceData
-from ..numeric_utils import _fast_kde, _fast_kde_2d
+from ..numeric_utils import _fast_kde_2d
+from ..kde_utils import kde
 from .plot_utils import get_plotting_function
 from ..rcparams import rcParams
 
@@ -14,7 +15,7 @@ def plot_kde(
     cumulative=False,
     rug=False,
     label=None,
-    bw=4.5,
+    bw=1,
     quantiles=None,
     rotated=False,
     contour=True,
@@ -50,17 +51,18 @@ def plot_kde(
         If True adds a rugplot. Defaults to False. Ignored for 2D KDE
     label : string
         Text to include as part of the legend
-    bw : float
-        Bandwidth scaling factor for 1D KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's
-        rule of thumb (the default rule used by SciPy).
+    bw: float, optional
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
     quantiles : list
-        Quantiles in ascending order used to segment the KDE. Use [.25, .5, .75] for quartiles.
-        Defaults to None.
+        Quantiles in ascending order used to segment the KDE.
+        Use [.25, .5, .75] for quartiles. Defaults to None.
     rotated : bool
         Whether to rotate the 1D KDE plot 90 degrees.
     contour : bool
-        If True plot the 2D KDE using contours, otherwise plot a smooth 2D KDE. Defaults to True.
+        If True plot the 2D KDE using contours, otherwise plot a smooth 2D KDE.
+        Defaults to True.
     fill_last : bool
         If True fill the last contour of the 2D KDE plot. Defaults to False.
     textsize: float
@@ -192,7 +194,8 @@ def plot_kde(
         )
 
     if values2 is None:
-        density, lower, upper = _fast_kde(values, cumulative, bw)
+        grid, density = kde(values, bw_fct=bw, cumulative=cumulative)
+        lower, upper = grid[0], grid[-1]
 
         if cumulative:
             density_q = density
diff --git a/arviz/plots/loopitplot.py b/arviz/plots/loopitplot.py
index 9d6268fa32..9d4823a042 100644
--- a/arviz/plots/loopitplot.py
+++ b/arviz/plots/loopitplot.py
@@ -3,7 +3,7 @@
 import scipy.stats as stats
 
 from ..stats import loo_pit as _loo_pit
-from ..numeric_utils import _fast_kde
+from ..kde_utils import kde
 from .plot_utils import get_plotting_function
 
 from ..rcparams import rcParams
@@ -161,10 +161,9 @@ def plot_loo_pit(
         )
         unif_ecdf = unif_ecdf / n_data_points
     else:
-        loo_pit_kde, xmin, xmax = _fast_kde(loo_pit)
+        x_vals, loo_pit_kde = kde(loo_pit)
 
         unif = np.random.uniform(size=(n_unif, loo_pit.size))
-        x_vals = np.linspace(xmin, xmax, len(loo_pit_kde))
         if use_hdi:
             n_obs = loo_pit.size
             hdi_ = stats.beta(n_obs / 2, n_obs / 2).ppf((1 - credible_interval) / 2)
diff --git a/arviz/plots/plot_utils.py b/arviz/plots/plot_utils.py
index 2a52cab39c..6aad9984f8 100644
--- a/arviz/plots/plot_utils.py
+++ b/arviz/plots/plot_utils.py
@@ -13,7 +13,7 @@
 from scipy.stats import mode
 import xarray as xr
 
-from ..numeric_utils import _fast_kde
+from ..kde_utils import kde
 from ..rcparams import rcParams
 
 KwargSpec = Dict[str, Any]
@@ -611,7 +611,7 @@ def get_plotting_function(plot_name, plot_module, backend):
     return plotting_method
 
 
-def calculate_point_estimate(point_estimate, values, bw=4.5):
+def calculate_point_estimate(point_estimate, values, bw=1):
     """Validate and calculate the point estimate.
 
     Parameters
@@ -620,10 +620,10 @@ def calculate_point_estimate(point_estimate, values, bw=4.5):
         Plot point estimate per variable. Values should be 'mean', 'median', 'mode' or None.
         Defaults to 'auto' i.e. it falls back to default set in rcParams.
     values : 1-d array
-    bw : float
-        Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the
-        KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb
-        (the default used rule by SciPy).
+    bw: float, optional
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
 
     Returns
     -------
@@ -643,8 +643,7 @@ def calculate_point_estimate(point_estimate, values, bw=4.5):
         point_value = values.mean()
     elif point_estimate == "mode":
         if isinstance(values[0], float):
-            density, lower, upper = _fast_kde(values, bw=bw)
-            x = np.linspace(lower, upper, len(density))
+            x, density = kde(values, bw_fct=bw)
             point_value = x[np.argmax(density)]
         else:
             point_value = mode(values)[0][0]
@@ -713,8 +712,7 @@ def sample_reference_distribution(dist, shape):
     densities = []
     dist_rvs = dist.rvs(size=shape)
     for idx in range(shape[1]):
-        density, xmin, xmax = _fast_kde(dist_rvs[:, idx])
-        x_s = np.linspace(xmin, xmax, len(density))
+        x_s, density = kde(dist_rvs[:, idx])
         x_ss.append(x_s)
         densities.append(density)
     return np.array(x_ss).T, np.array(densities).T
diff --git a/arviz/plots/posteriorplot.py b/arviz/plots/posteriorplot.py
index 04d0307f6a..d502c07d34 100644
--- a/arviz/plots/posteriorplot.py
+++ b/arviz/plots/posteriorplot.py
@@ -26,7 +26,7 @@ def plot_posterior(
     rope=None,
     ref_val=None,
     kind="kde",
-    bw=4.5,
+    bw=1,
     bins=None,
     ax=None,
     backend=None,
@@ -82,10 +82,11 @@ def plot_posterior(
     kind: str
         Type of plot to display (kde or hist) For discrete variables this argument is ignored and
         a histogram is always used.
-    bw: float
-        Bandwidth scaling factor for the KDE. Should be larger than 0. The higher this number the
-        smoother the KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule
-        of thumb (the default rule used by SciPy). Only works if `kind == kde`.
+    bw : float
+        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
+        The higher this number the smoother the KDE will be.
+        Defaults to 1 which means the bandwidth is not modified.
+        Only works if `kind == kde`.
     bins: integer or sequence or 'auto', optional
         Controls the number of bins, accepts the same keywords `matplotlib.hist()` does. Only works
         if `kind == hist`. If None (default) it will use `auto` for continuous variables and
diff --git a/arviz/plots/violinplot.py b/arviz/plots/violinplot.py
index e69dd9012b..81fa9a868c 100644
--- a/arviz/plots/violinplot.py
+++ b/arviz/plots/violinplot.py
@@ -19,7 +19,7 @@ def plot_violin(
     rug=False,
     hdi_prob=None,
     shade=0.35,
-    bw=4.5,
+    bw=1,
     sharex=True,
     sharey=True,
     figsize=None,
@@ -64,9 +64,8 @@ def plot_violin(
         Alpha blending value for the shaded area under the curve, between 0
         (no shade) and 1 (opaque). Defaults to 0
     bw: float
-        Bandwidth scaling factor. Should be larger than 0. The higher this number the smoother the
-        KDE will be. Defaults to 4.5 which is essentially the same as the Scott's rule of thumb
-        (the default rule used by SciPy).
+        Bandwidth scaling factor. Must be larger than 0. The higher this number the smoother the
+        KDE will be. Defaults to 1 (no modification)
     figsize: tuple
         Figure size. If None it will be defined automatically.
     textsize: int
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index f7dd0b4639..94a7a1072d 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -23,7 +23,8 @@
     _circular_standard_deviation,
     get_log_likelihood as _get_log_likelihood,
 )
-from ..numeric_utils import _fast_kde, histogram, get_bins
+from ..numeric_utils import histogram, get_bins
+from ..kde_utils import kde
 from ..utils import _var_names, Numba, _numba_var, get_coords, credible_interval_warning
 from ..rcparams import rcParams
 
@@ -545,10 +546,10 @@ def _hdi_multimodal(ary, hdi_prob, skipna, max_modes):
         ary = ary[~np.isnan(ary)]
 
     if ary.dtype.kind == "f":
-        density, lower, upper = _fast_kde(ary)
+        bins, density = kde(ary)
+        lower, upper = bins[0], bins[-1]
         range_x = upper - lower
         dx = range_x / len(density)
-        bins = np.linspace(lower, upper, len(density))
     else:
         bins = get_bins(ary)
         _, density, _ = histogram(ary, bins=bins)
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index b6e6d49064..b7209e4de7 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -45,7 +45,7 @@
     plot_bpv,
 )
 from ...utils import _cov
-from ...numeric_utils import _fast_kde
+from ...kde_utils import kde
 
 rcParams["data.load"] = "eager"
 
@@ -923,24 +923,25 @@ def test_plot_hdi_dataset_error(models):
 
 
 @pytest.mark.parametrize("limits", [(-10.0, 10.0), (-5, 5), (None, None)])
-def test_fast_kde_scipy(limits):
+def test_kde_scipy(limits):
+    """
+    Evaluates if sum of density is the same for our implementation
+    and the implementation in scipy
+    """
     data = np.random.normal(0, 1, 10000)
-    if limits[0] is None:
-        x = np.linspace(data.min(), data.max(), 200)  # pylint: disable=no-member
-    else:
-        x = np.linspace(*limits, 500)
-    density = gaussian_kde(data).evaluate(x)
-    density_fast = _fast_kde(data, xmin=limits[0], xmax=limits[1])[0]
-
-    np.testing.assert_almost_equal(density_fast.sum(), density.sum(), 1)
+    grid, density_own = kde(data, custom_lims=limits)[1]
+    density_sp = gaussian_kde(data).evaluate(grid)
+    np.testing.assert_almost_equal(density_own.sum(), density_sp.sum(), 1)
 
 
 @pytest.mark.parametrize("limits", [(-10.0, 10.0), (-5, 5), (None, None)])
-def test_fast_kde_cumulative(limits):
+def test_kde_cumulative(limits):
+    """
+    Evaluates if last value of cumulative density is 1
+    """
     data = np.random.normal(0, 1, 1000)
-    density_fast = _fast_kde(data, xmin=limits[0], xmax=limits[1], cumulative=True)[0]
-    np.testing.assert_almost_equal(round(density_fast[-1], 3), 1)
-
+    density = kde(data, custom_lims=limits, cumulative=True)[1]
+    np.testing.assert_almost_equal(round(density[-1], 3), 1)
 
 @pytest.mark.parametrize(
     "kwargs",

From 749976429c45a465bc10536bafea73b37faddcc1 Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Fri, 17 Jul 2020 20:58:46 -0300
Subject: [PATCH 02/11] Update according to requests

---
 arviz/kde_utils.py | 37 +++++++++++++++++--------------------
 1 file changed, 17 insertions(+), 20 deletions(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index f83a6bc07f..357eadc2f1 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -1,14 +1,13 @@
 """Kernel density estimation functions for ArviZ."""
 
+import warnings
 import numpy as np
-from scipy.signal import gaussian, convolve
 from scipy.fftpack import fft
 from scipy.optimize import fsolve
 from scipy.optimize import brentq
+from scipy.signal import gaussian, convolve
 from scipy.special import ive
-import warnings
 
-# Own
 from .stats.stats_utils import histogram
 
 # Bandwidth functions ---------------------------------------------------------------
@@ -40,8 +39,9 @@ def bw_isj(x, grid_counts=None, x_range=None):
     """
     Improved Sheather and Jones method as explained in [1]
     This is an internal version pretended to be used by the KDE estimator.
-    When used internally computation time is saved because
-    things like minimums, maximums and the grid are pre-computed.
+    When used internally computation time is saved because things like minimums,
+
+    maximums and the grid are pre-computed.
 
     References
     ----------
@@ -85,7 +85,7 @@ def bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
     return 0.5 * (bw_silverman(x, x_std) + bw_isj(x, grid_counts, x_range))
 
 
-def bw_circular(x):
+def bw_taylor(x):
     """
     Bandwidth selector for circular kernel density estimation
     as introduced in [1].
@@ -205,12 +205,12 @@ def circular_mean(x):
     return mean
 
 
-def A1inv(x):
+def a1inv(x):
     """
     Inverse function of the ratio of the first and zeroth order
     Bessel functions of the first kind.
 
-    Returns the value k, such that A1inv(x) = k, i.e. A1(k) = x.
+    Returns the value k, such that a1inv(x) = k, i.e. a1(k) = x.
     """
     if 0 <= x < 0.53:
         return 2 * x + x ** 3 + (5 * x ** 5) / 6
@@ -222,7 +222,7 @@ def A1inv(x):
 
 def kappa_mle(x):
     mean = circular_mean(x)
-    kappa = A1inv(np.mean(np.cos(x - mean)))
+    kappa = a1inv(np.mean(np.cos(x - mean)))
     return kappa
 
 
@@ -397,8 +397,7 @@ def check_custom_lims(custom_lims, x_min, x_max):
     if custom_lims[1] is None:
         custom_lims[1] = x_max
 
-    types = (int, float, np.integer, np.float)
-    all_numeric = all(isinstance(i, types) for i in custom_lims)
+    all_numeric = all(isinstance(i, (int, float, np.integer, np.float)) for i in custom_lims)
     if not all_numeric:
         raise TypeError((
             f"Elements of `custom_lims` must be numeric or None.\n"
@@ -452,9 +451,7 @@ def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True
 
     """
 
-    # Set up number of bins. Should enable larger grids?
-    if grid_len > 4096:
-        grid_len = 4096
+    # Set up number of bins.
     if grid_len < 100:
         grid_len = 100
     grid_len = int(grid_len)
@@ -619,7 +616,7 @@ def kde_linear(
 
 def kde_circular(
     x,
-    bw="default",
+    bw="taylor",
     bw_fct=1,
     bw_return=False,
     custom_lims=None,
@@ -642,7 +639,7 @@ def kde_circular(
     bw: int, float or str, optional
         If numeric, indicates the bandwidth and must be positive.
         If str, indicates the method to estimate the bandwidth and must be
-        one "default". Defaults to "default".
+        "taylor" since it is the only option supported so far. Defaults to "taylor".
     bw_fct: float, optional
         A value that multiplies `bw` which enables tuning smoothness by hand.
         Must be positive. Values above 1 decrease smoothness while values
@@ -676,17 +673,17 @@ def kde_circular(
     if isinstance(bw, bool):
         raise ValueError((
             "`bw` can't be of type `bool`.\n"
-            "Expected a positive numeric or 'default'"
+            "Expected a positive numeric or 'taylor'"
         ))
     if isinstance(bw, (int, float)):
         if bw < 0:
             raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
     if isinstance(bw, str):
-        if bw == "default":
-            bw = bw_circular(x)
+        if bw == "taylor":
+            bw = bw_taylor(x)
         else:
             raise ValueError((
-                f"`bw` must be a positive numeric or `default`, not {bw}"
+                f"`bw` must be a positive numeric or `taylor`, not {bw}"
             ))
     bw *= bw_fct
 

From 8335720da23327abb8f640703b21419c7dfc09ad Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Thu, 23 Jul 2020 01:07:44 -0300
Subject: [PATCH 03/11] Fix how circular data is handled and listed changes in
 CHANGELOG. * In the previous version the function  tried to normalize the
 angles to [-pi, pi) but it was not correct. This commit includes a function 
 that works appropiately.

---
 CHANGELOG.md       |   7 +++
 arviz/kde_utils.py | 106 ++++++++++++++++++++++++++++++++++++++++-----
 2 files changed, 103 insertions(+), 10 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b62a3cca88..1f0cd5808c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,13 @@
 * update `radon` example dataset to current InferenceData schema specification (#1320)
 * update `from_cmdstan` functionality and add warmup groups (#1330)
 * restructure plotting code to be compatible with mpl>=3.3 (#1312)
+* Replaced `_fast_kde()` with `kde()` which now also supports circular variables
+via the argument `circular` (#1284).
+
+### Maintenance and fixes
+* plot_posterior: fix overlap of hdi and rope (#1263)
+* All the functions that used to call `_fast_kde`() now use `kde()` and have
+been updated to handle the new types returned.
 
 ### Deprecation
 
diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 357eadc2f1..7033d836ce 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -40,7 +40,6 @@ def bw_isj(x, grid_counts=None, x_range=None):
     Improved Sheather and Jones method as explained in [1]
     This is an internal version pretended to be used by the KDE estimator.
     When used internally computation time is saved because things like minimums,
-
     maximums and the grid are pre-computed.
 
     References
@@ -111,7 +110,7 @@ def bw_taylor(x):
     return (num / den) ** 0.4
 
 
-BW_METHODS = {
+BW_METHODS_LINEAR = {
         "scott": bw_scott,
         "silverman": bw_silverman,
         "isj": bw_isj,
@@ -137,13 +136,13 @@ def select_bw_method(method):
     """
 
     method_lower = method.lower()
-    if method_lower not in BW_METHODS.keys():
+    if method_lower not in BW_METHODS_LINEAR.keys():
         raise ValueError((
             f"Unrecognized bandwidth method.\n"
             f"Input is: {method}.\n"
-            f"Expected one of: {list(BW_METHODS.keys())}."
+            f"Expected one of: {list(BW_METHODS_LINEAR.keys())}."
         ))
-    bw_fun = BW_METHODS[method_lower]
+    bw_fun = BW_METHODS_LINEAR[method_lower]
     return bw_fun
 
 
@@ -170,7 +169,7 @@ def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
         raise ValueError((
             f"`bw` must not be of type `bool`.\n"
             f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(BW_METHODS.keys())}."))
+            f"{list(BW_METHODS_LINEAR.keys())}."))
     if isinstance(bw, (int, float)):
         if bw < 0:
             raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
@@ -187,11 +186,21 @@ def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
         raise ValueError((
             f"Unrecognized `bw` argument.\n"
             f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(BW_METHODS.keys())}."))
+            f"{list(BW_METHODS_LINEAR.keys())}."))
     return bw
 
 
-# Misc utils
+# Misc utils ------------------------------------------------------------------------
+def normalize_angle(x, pi_centered=False):
+    """
+    Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
+    depending on `pi_centered`.
+    """
+    if pi_centered:
+        return x % (2 * np.pi)
+    else:
+        return (x + np.pi) % (2 * np.pi) - np.pi
+    
 def vonmises_pdf(x, mu, kappa):
     assert kappa > 0, "Argument 'kappa' must be positive."
     pdf = 1 / (2 * np.pi * ive(0, kappa)) * np.exp(np.cos(x - mu) - 1) ** kappa
@@ -495,6 +504,84 @@ def kde(x, circular=False, **kwargs):
     grid : Gridded numpy array for the x values.
     pdf : Numpy array for the density estimates.
     bw: optional, the estimated bandwidth.
+    
+    Examples
+    --------
+    Default density estimation for linear data
+    .. plot::
+        :context: close-figs
+        
+    	>>> import numpy as np
+	>>> import matplotlib.pyplot as plt
+	>>> from arviz.kde_utils import kde
+	>>> 
+	>>> rvs = np.random.gamma(shape=1.8, size=1000)
+	>>> grid, pdf = kde(rvs)
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Density estimation for linear data with Silverman's rule bandwidth
+    .. plot::
+        :context: close-figs
+
+	>>> grid, pdf = kde(rvs, bw="silverman")
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Density estimation for linear data with scaled bandwidth
+    .. plot::
+        :context: close-figs
+
+	>>> # bw_fct > 1 means more smoothness.
+	>>> grid, pdf = kde(rvs, bw_fct=2.5)
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Default density estimation for linear data with extended limits
+    .. plot::
+        :context: close-figs
+
+	>>> grid, pdf = kde(rvs, bound_correction=False, extend=True, extend_fct=0.5)
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Default density estimation for linear data with custom limits
+    .. plot::
+        :context: close-figs
+	# It accepts tuples and lists of length 2.
+	>>> grid, pdf = kde(rvs, bound_correction=False, custom_lims=(0, 10))
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Default density estimation for circular data
+    .. plot::
+        :context: close-figs
+
+	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
+	>>> grid, pdf = kde(rvs, circular=True)
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Density estimation for circular data with scaled bandwidth
+    .. plot::
+        :context: close-figs
+
+	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
+	>>> # bw_fct > 1 means less smoothness.
+	>>> grid, pdf = kde(rvs, circular=True, bw_fct=3)
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+    Density estimation for circular data with custom limits
+    .. plot::
+        :context: close-figs
+	>>> # This is still experimental, does not always work.
+	>>> rvs = np.random.vonmises(mu=0, kappa=30, size=500)
+	>>> grid, pdf = kde(rvs, circular=True, custom_lims=(-1, 1))
+	>>> plt.plot(grid, pdf)
+	>>> plt.show()
+	
+	
     """
     if circular:
         kde_fun = kde_circular
@@ -662,8 +749,7 @@ def kde_circular(
     x = check_type(x)
 
     # All values between -pi and pi
-    x[x > np.pi] = x[x > np.pi] - 2 * np.pi
-    x[x < -np.pi] = x[x < -np.pi] + 2 * np.pi
+    x = normalize_angle(x)
 
     # Assert `bw_fct` is numeric and positive
     assert isinstance(bw_fct, (int, float))

From eaaff9d045fe1d71f7538db0fd462ffc7b75122d Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Thu, 23 Jul 2020 10:23:43 -0300
Subject: [PATCH 04/11] Update according to required changes.

---
 CHANGELOG.md       | 4 ++--
 arviz/kde_utils.py | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f0cd5808c..d07ef7752c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,11 +15,11 @@
 * restructure plotting code to be compatible with mpl>=3.3 (#1312)
 * Replaced `_fast_kde()` with `kde()` which now also supports circular variables
 via the argument `circular` (#1284).
+* Replaced `_fast_kde()` with `kde()` which now also supports circular variables via the argument `circular` (#1284).
 
 ### Maintenance and fixes
 * plot_posterior: fix overlap of hdi and rope (#1263)
-* All the functions that used to call `_fast_kde`() now use `kde()` and have
-been updated to handle the new types returned.
+* All the functions that used to call `_fast_kde`() now use `kde()` and have been updated to handle the new types returned.
 
 ### Deprecation
 
diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 7033d836ce..f700055e49 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -10,7 +10,7 @@
 
 from .stats.stats_utils import histogram
 
-# Bandwidth functions ---------------------------------------------------------------
+# Bandwidth functions
 # Linear KDE
 def bw_scott(x, x_std=None):
     """
@@ -190,7 +190,7 @@ def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
     return bw
 
 
-# Misc utils ------------------------------------------------------------------------
+# Misc utils
 def normalize_angle(x, pi_centered=False):
     """
     Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
@@ -330,7 +330,7 @@ def _root(function, N, args, x):
             bw = (bw_silverman(x) / np.ptp(x)) ** 2
     return bw
 
-# KDE Utilities ---------------------------------------------------------------------
+# KDE Utilities
 
 def check_type(x):
     """
@@ -481,7 +481,7 @@ def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True
         grid_max = x_max
     return grid_min, grid_max, grid_len
 
-# KDE Functions --------------------------------------------------------------------
+# KDE Functions
 
 def kde(x, circular=False, **kwargs):
     """

From 5838824c394d7b3d80bf2e65051483b85e9a4f21 Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Wed, 29 Jul 2020 19:58:11 -0300
Subject: [PATCH 05/11] Several changes according to the PR discussion and
 other suggestions received. * Functions _circular_mean() and
 normalize_angle() now are in numeric_utils.py. * If the array has more than 1
 dimension, _check_type() tries to flatten it silenty. If it is not possible,
 an error is raised. * Improved how bw_fct type and value are checked and
 reported. * All functions now start with _ because they are all internals. *
 _select_bw_method() no longer exists. Its functionality is now incorporated
 into _get_bw(). * Changed internal _bw_* functions to a common signature. *
 Added DeprecationWarning to _fast_kde() and now it calls _kde(). * Added
 tests for _normalize_angle() and _circular_mean(). The following functions
 and their backends have been modified to incorporate the circular argument
 that is passed to _kde() as well as the new meaning of the argument bw (see
 below): * plot_density() * plot_kde() * plot_dist() * plot_posterior() *
 plot_violion() * calculate_point_estimate() In addition, plot_kde() gained a
 new argument, adaptive, which is used to determine if the bandwidth is
 adaptive or not. Finally, note the meaning of the argument bw has changed. It
 used to be a bandwidth multiplication factor Now it can be either a string
 indicating the bandwidth method or the bandwidth itself. The internal
 function _kde() has an argument bw_fct which is the multiplicative factor.

---
 CHANGELOG.md                                  |   2 +-
 arviz/kde_utils.py                            | 335 ++++++++----------
 arviz/numeric_utils.py                        |  41 ++-
 arviz/plots/backends/bokeh/bpvplot.py         |   6 +-
 arviz/plots/backends/bokeh/densityplot.py     |   9 +-
 arviz/plots/backends/bokeh/distplot.py        |   2 +
 arviz/plots/backends/bokeh/forestplot.py      |   4 +-
 arviz/plots/backends/bokeh/loopitplot.py      |   4 +-
 arviz/plots/backends/bokeh/posteriorplot.py   |   6 +-
 arviz/plots/backends/bokeh/ppcplot.py         |   4 +-
 arviz/plots/backends/bokeh/violinplot.py      |  16 +-
 arviz/plots/backends/matplotlib/bpvplot.py    |   6 +-
 .../plots/backends/matplotlib/densityplot.py  |  22 +-
 arviz/plots/backends/matplotlib/distplot.py   |   2 +
 arviz/plots/backends/matplotlib/forestplot.py |   4 +-
 arviz/plots/backends/matplotlib/loopitplot.py |   5 +-
 .../backends/matplotlib/posteriorplot.py      |   6 +-
 arviz/plots/backends/matplotlib/ppcplot.py    |   6 +-
 arviz/plots/backends/matplotlib/violinplot.py |  14 +-
 arviz/plots/densityplot.py                    |  24 +-
 arviz/plots/distplot.py                       |  18 +-
 arviz/plots/energyplot.py                     |  10 +-
 arviz/plots/kdeplot.py                        |  65 +++-
 arviz/plots/loopitplot.py                     |   4 +-
 arviz/plots/plot_utils.py                     |  27 +-
 arviz/plots/posteriorplot.py                  |  18 +-
 arviz/plots/violinplot.py                     |  17 +-
 arviz/stats/stats.py                          |   4 +-
 .../tests/base_tests/test_plots_matplotlib.py |   6 +-
 arviz/tests/base_tests/test_utils.py          |  22 ++
 30 files changed, 429 insertions(+), 280 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index d07ef7752c..b36ab8e1c3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,7 +19,7 @@ via the argument `circular` (#1284).
 
 ### Maintenance and fixes
 * plot_posterior: fix overlap of hdi and rope (#1263)
-* All the functions that used to call `_fast_kde`() now use `kde()` and have been updated to handle the new types returned.
+* All the functions that used to call `_fast_kde`() now use `kde()` and have been updated to handle the new types returned (#1284).
 
 ### Deprecation
 
diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index f700055e49..9ab599f64a 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -9,12 +9,11 @@
 from scipy.special import ive
 
 from .stats.stats_utils import histogram
+import arviz.numeric_utils as numeric_utils
+
+def _bw_scott(x, grid_counts=None, x_std=None, x_range=None):
+    """Scott's Rule.
 
-# Bandwidth functions
-# Linear KDE
-def bw_scott(x, x_std=None):
-    """
-    Scott's Rule
     """
     if x_std is None:
         x_std = np.std(x)
@@ -22,9 +21,9 @@ def bw_scott(x, x_std=None):
     return bw
 
 
-def bw_silverman(x, x_std=None):
-    """
-    Silverman's Rule.
+def _bw_silverman(x, grid_counts=None, x_std=None, x_range=None):
+    """Silverman's Rule.
+
     """
     if x_std is None:
         x_std = np.std(x)
@@ -35,9 +34,11 @@ def bw_silverman(x, x_std=None):
     return bw
 
 
-def bw_isj(x, grid_counts=None, x_range=None):
+def _bw_isj(x, grid_counts=None, x_std=None, x_range=None):
     """
-    Improved Sheather and Jones method as explained in [1]
+    Improved Sheather-Jones bandwidth estimation.
+
+    Improved Sheather and Jones method as explained in [1]_.
     This is an internal version pretended to be used by the KDE estimator.
     When used internally computation time is saved because things like minimums,
     maximums and the grid are pre-computed.
@@ -68,7 +69,7 @@ def bw_isj(x, grid_counts=None, x_range=None):
     grid_relfreq = grid_counts / x_len
 
     # Discrete cosine transform of the data
-    a_k = dct1d(grid_relfreq)
+    a_k = _dct1d(grid_relfreq)
 
     k_sq = np.arange(1, grid_len) ** 2
     a_sq = a_k[range(1, grid_len)] ** 2
@@ -77,20 +78,21 @@ def bw_isj(x, grid_counts=None, x_range=None):
     h = t ** 0.5 * x_range
     return h
 
-def bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
-    """
-    Experimental bandwidth estimator.
+def _bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
+    """Experimental bandwidth estimator.
+
     """
-    return 0.5 * (bw_silverman(x, x_std) + bw_isj(x, grid_counts, x_range))
+    bw_silverman = _bw_silverman(x, x_std=x_std)
+    bw_isj = _bw_isj(x, grid_counts=grid_counts, x_range=x_range)
+    return 0.5 * (bw_silverman + bw_isj)
 
 
-def bw_taylor(x):
-    """
-    Bandwidth selector for circular kernel density estimation
-    as introduced in [1].
+def _bw_taylor(x):
+    """Taylor's rule for circular bandwidth estimation.
+
     This function implements a rule-of-thumb for choosing the bandwidth of
     a von Mises kernel density estimator that assumes the underlying
-    distribution is von Mises.
+    distribution is von Mises as introduced in [1]_.
     It is analogous to Scott's rule for the Gaussian KDE.
 
     Circular bandwidth has a different scale from linear bandwidth.
@@ -99,54 +101,25 @@ def bw_taylor(x):
 
     References
     ----------
-    [1] C.C Taylor (2008). Automatic bandwidth selection
-    for circular density estimation.
-    Computational Statistics and Data Analysis, 52, 7, 3493–3500.
+    .. [1] C.C Taylor (2008). Automatic bandwidth selection for circular
+           density estimation.
+           Computational Statistics and Data Analysis, 52, 7, 3493–3500.
     """
     x_len = len(x)
-    kappa = kappa_mle(x)
+    kappa = _kappa_mle(x)
     num = 3 * x_len * kappa ** 2 * ive(2, 2 * kappa)
     den = 4 * np.pi ** 0.5 * ive(0, kappa) ** 2
     return (num / den) ** 0.4
 
 
-BW_METHODS_LINEAR = {
-        "scott": bw_scott,
-        "silverman": bw_silverman,
-        "isj": bw_isj,
-        "experimental": bw_experimental
+_BW_METHODS_LINEAR = {
+        "scott": _bw_scott,
+        "silverman": _bw_silverman,
+        "isj": _bw_isj,
+        "experimental": _bw_experimental
     }
 
-
-def select_bw_method(method):
-    """
-    Selects a function to compute the bandwidth.
-    Also checks method `bw` is correctly specified.
-    Otherwise, throws an error.
-
-    Parameters
-    ----------
-    method : str
-        Method to estimate the bandwidth.
-
-    Returns
-    -------
-    bw_fun: function
-        Function to compute the bandwidth.
-    """
-
-    method_lower = method.lower()
-    if method_lower not in BW_METHODS_LINEAR.keys():
-        raise ValueError((
-            f"Unrecognized bandwidth method.\n"
-            f"Input is: {method}.\n"
-            f"Expected one of: {list(BW_METHODS_LINEAR.keys())}."
-        ))
-    bw_fun = BW_METHODS_LINEAR[method_lower]
-    return bw_fun
-
-
-def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
+def _get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
     """
     Computes bandwidth for a given data `x` and `bw`.
     Also checks `bw` is correctly specified.
@@ -169,55 +142,41 @@ def get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
         raise ValueError((
             f"`bw` must not be of type `bool`.\n"
             f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(BW_METHODS_LINEAR.keys())}."))
+            f"{list(_BW_METHODS_LINEAR.keys())}."))
     if isinstance(bw, (int, float)):
         if bw < 0:
             raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
     elif isinstance(bw, str):
-        bw_fun = select_bw_method(bw)
         bw_lower = bw.lower()
-        if bw_lower == "isj":
-            bw = bw_fun(x, grid_counts, x_range)
-        elif bw_lower in ["scott", "silverman"]:
-            bw = bw_fun(x, x_std)
-        elif bw_lower == "experimental":
-            bw = bw_fun(x, grid_counts, x_std, x_range)
+
+        if bw_lower not in _BW_METHODS_LINEAR.keys():
+            raise ValueError((
+                f"Unrecognized bandwidth method.\n"
+                f"Input is: {method}.\n"
+                f"Expected one of: {list(_BW_METHODS_LINEAR.keys())}."
+            ))
+
+        bw_fun = _BW_METHODS_LINEAR[bw_lower]
+        bw = bw_fun(x, grid_counts, x_std, x_range)           
     else:
         raise ValueError((
             f"Unrecognized `bw` argument.\n"
             f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(BW_METHODS_LINEAR.keys())}."))
+            f"{list(_BW_METHODS_LINEAR.keys())}."))
     return bw
 
 
-# Misc utils
-def normalize_angle(x, pi_centered=False):
-    """
-    Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
-    depending on `pi_centered`.
-    """
-    if pi_centered:
-        return x % (2 * np.pi)
-    else:
-        return (x + np.pi) % (2 * np.pi) - np.pi
-    
-def vonmises_pdf(x, mu, kappa):
-    assert kappa > 0, "Argument 'kappa' must be positive."
+def _vonmises_pdf(x, mu, kappa):
+    if kappa <= 0:
+        raise ValueError("Argument 'kappa' must be positive.")
     pdf = 1 / (2 * np.pi * ive(0, kappa)) * np.exp(np.cos(x - mu) - 1) ** kappa
     return pdf
 
 
-def circular_mean(x):
-    sinr = np.sum(np.sin(x))
-    cosr = np.sum(np.cos(x))
-    mean = np.arctan2(sinr, cosr)
-    return mean
-
-
-def a1inv(x):
+def _a1inv(x):
     """
-    Inverse function of the ratio of the first and zeroth order
-    Bessel functions of the first kind.
+    Inverse function of the ratio of the first and
+    zeroth order Bessel functions of the first kind.
 
     Returns the value k, such that a1inv(x) = k, i.e. a1(k) = x.
     """
@@ -229,13 +188,13 @@ def a1inv(x):
         return 1 / (x ** 3 - 4 * x ** 2 + 3 * x)
 
 
-def kappa_mle(x):
-    mean = circular_mean(x)
-    kappa = a1inv(np.mean(np.cos(x - mean)))
+def _kappa_mle(x):
+    mean = numeric_utils._circular_mean(x)
+    kappa = _a1inv(np.mean(np.cos(x - mean)))
     return kappa
 
 
-def dct1d(x):
+def _dct1d(x):
     """
     Discrete Cosine Transform in 1 Dimension
 
@@ -268,8 +227,7 @@ def dct1d(x):
 
 def _fixed_point(t, N, k_sq, a_sq):
     """
-    Implementation of the function t-zeta*gamma^[l](t) derived
-    from equation (30) in [1]
+    Implementation of the function t-zeta*gamma^[l](t) derived from equation (30) in [1]
 
     References
     ----------
@@ -298,45 +256,35 @@ def _fixed_point(t, N, k_sq, a_sq):
 
 
 def _root(function, N, args, x):
-    # The idea here was borrowed from KDEpy implementation.
     # The right bound is at most 0.01
-    found = 0
+    found = False
     N = max(min(1050, N), 50)
     tol = 10e-12 + 0.01 * (N - 50) / 1000
 
-    while found == 0:
+    while not found:
         try:
-            bw, res = brentq(function, 0, 0.01, args=args, full_output=True,
-                             disp=False)
-            found = 1 if res.converged else 0
+            bw, res = brentq(function, 0, 0.01, args=args, full_output=True, disp=False)
+            found = True if res.converged else False
         except ValueError:
             bw = 0
             tol *= 2.0
-            found = 0
-        if bw <= 0:
-            warnings.warn(
-                "Improved Sheather-Jones did not converge to a positive value. "
-                "Using Silverman's rule instead.",
-                Warning
-            )
-            bw = (bw_silverman(x) / np.ptp(x)) ** 2
+            found = False
+        if bw <= 0 or tol >= 1:
+            # warnings.warn(
+            #    "Improved Sheather-Jones did not converge as expected. "
+            #    "Using Silverman's rule instead.",
+            #    Warning
+            # )
+            bw = (_bw_silverman(x) / np.ptp(x)) ** 2
             return bw
-        if tol >= 1:
-            warnings.warn(
-                "Improved Sheather-Jones did not converge. "
-                "Using Silverman's rule instead.",
-                Warning
-            )
-            bw = (bw_silverman(x) / np.ptp(x)) ** 2
     return bw
 
-# KDE Utilities
 
-def check_type(x):
+def _check_type(x):
     """
     Checks the input is of the correct type.
-    It only accepts numeric lists/numpy arrays of 1 dimension.
-    If input is not of the correct type, an informative message is thrown.
+    It only accepts numeric lists/numpy arrays of 1 dimension or something that
+    can be flattened to 1 dimension.
 
     Parameters
     ----------
@@ -351,27 +299,27 @@ def check_type(x):
     """
     if not isinstance(x, (list, np.ndarray)):
         raise ValueError((
-            f"`x` is of the wrong type.\n"
             f"Can't produce a density estimator for {type(x)}.\n"
-            f"Please input a numeric list or numpy array."
+            f"Please input a list with numbers or numpy array."
         ))
-
-    # Will raise an error if `x` can't be casted to numeric
-    x = np.asfarray(x)
+    # Will raise an error if `x` can't be casted to numeric or flattened to one dimension.
+    try:
+        x = np.asfarray(x).flatten()
+    except Exception as e:
+        print(
+            "The following exception occurred while trying to convert `x`",
+            "to a 1 dimensional float array."
+        )
+        raise e
 
     x = x[np.isfinite(x)]
     if x.size == 0:
         raise ValueError("`x` does not contain any finite number.")
 
-    if x.ndim != 1:
-        raise ValueError((
-            f"Unsupported dimension number.\n"
-            f"Density estimator only works with 1-dimensional data, not {x.ndim}."
-        ))
     return x
 
 
-def check_custom_lims(custom_lims, x_min, x_max):
+def _check_custom_lims(custom_lims, x_min, x_max):
     """
     Checks whether `custom_lims` are of the correct type.
     It accepts numeric lists/tuples of length 2.
@@ -419,8 +367,8 @@ def check_custom_lims(custom_lims, x_min, x_max):
 
     return custom_lims
 
-def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True,
-            bound_correction=False):
+def _get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True,
+              bound_correction=False):
     """
     Computes the grid that bins the data used to estimate the density function
 
@@ -469,7 +417,7 @@ def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True
     # `custom_lims` overrides `extend`
     # `bound_correction` overrides `extend`
     if custom_lims is not None:
-        custom_lims = check_custom_lims(custom_lims, x_min, x_max)
+        custom_lims = _check_custom_lims(custom_lims, x_min, x_max)
         grid_min = custom_lims[0]
         grid_max = custom_lims[1]
     elif extend and not bound_correction:
@@ -481,11 +429,10 @@ def get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True
         grid_max = x_max
     return grid_min, grid_max, grid_len
 
-# KDE Functions
 
-def kde(x, circular=False, **kwargs):
-    """
-    1 dimensional density estimation.
+def _kde(x, circular=False, **kwargs):
+    """1 dimensional density estimation.
+
     It is a wrapper around `kde_linear()` and `kde_circular()`.
 
     Parameters
@@ -504,94 +451,98 @@ def kde(x, circular=False, **kwargs):
     grid : Gridded numpy array for the x values.
     pdf : Numpy array for the density estimates.
     bw: optional, the estimated bandwidth.
-    
+
     Examples
     --------
     Default density estimation for linear data
     .. plot::
         :context: close-figs
-        
+
     	>>> import numpy as np
 	>>> import matplotlib.pyplot as plt
-	>>> from arviz.kde_utils import kde
-	>>> 
+	>>> from arviz.kde_utils import _kde
+	>>>
 	>>> rvs = np.random.gamma(shape=1.8, size=1000)
-	>>> grid, pdf = kde(rvs)
+	>>> grid, pdf = _kde(rvs)
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Density estimation for linear data with Silverman's rule bandwidth
     .. plot::
         :context: close-figs
 
-	>>> grid, pdf = kde(rvs, bw="silverman")
+	>>> grid, pdf = _kde(rvs, bw="silverman")
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Density estimation for linear data with scaled bandwidth
     .. plot::
         :context: close-figs
 
 	>>> # bw_fct > 1 means more smoothness.
-	>>> grid, pdf = kde(rvs, bw_fct=2.5)
+	>>> grid, pdf = _kde(rvs, bw_fct=2.5)
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Default density estimation for linear data with extended limits
     .. plot::
         :context: close-figs
 
-	>>> grid, pdf = kde(rvs, bound_correction=False, extend=True, extend_fct=0.5)
+	>>> grid, pdf = _kde(rvs, bound_correction=False, extend=True, extend_fct=0.5)
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Default density estimation for linear data with custom limits
     .. plot::
         :context: close-figs
 	# It accepts tuples and lists of length 2.
-	>>> grid, pdf = kde(rvs, bound_correction=False, custom_lims=(0, 10))
+	>>> grid, pdf = __kde(rvs, bound_correction=False, custom_lims=(0, 10))
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Default density estimation for circular data
     .. plot::
         :context: close-figs
 
 	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
-	>>> grid, pdf = kde(rvs, circular=True)
+	>>> grid, pdf = _kde(rvs, circular=True)
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Density estimation for circular data with scaled bandwidth
     .. plot::
         :context: close-figs
 
 	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
 	>>> # bw_fct > 1 means less smoothness.
-	>>> grid, pdf = kde(rvs, circular=True, bw_fct=3)
+	>>> grid, pdf = _kde(rvs, circular=True, bw_fct=3)
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
+
     Density estimation for circular data with custom limits
     .. plot::
         :context: close-figs
 	>>> # This is still experimental, does not always work.
 	>>> rvs = np.random.vonmises(mu=0, kappa=30, size=500)
-	>>> grid, pdf = kde(rvs, circular=True, custom_lims=(-1, 1))
+	>>> grid, pdf = _kde(rvs, circular=True, custom_lims=(-1, 1))
 	>>> plt.plot(grid, pdf)
 	>>> plt.show()
-	
-	
+
+    See Also
+    --------
+    plot_kde : Compute and plot a kernel density estimate.
+    arviz.kde_utils._kde: Arviz KDE estimator
+
     """
     if circular:
-        kde_fun = kde_circular
+        kde_fun = _kde_circular
     else:
-        kde_fun = kde_linear
+        kde_fun = _kde_linear
 
     return kde_fun(x, **kwargs)
 
 
-def kde_linear(
+def _kde_linear(
     x,
     bw="experimental",
     adaptive=False,
@@ -602,7 +553,8 @@ def kde_linear(
     bw_return=False,
     custom_lims=None,
     cumulative=False,
-    grid_len=512
+    grid_len=512,
+    **kwargs
 ):
     """
     1 dimensional density estimation for linear data.
@@ -661,11 +613,14 @@ def kde_linear(
     """
 
     # Check `x` is from appropiate type
-    x = check_type(x)
+    x = _check_type(x)
+
+    # Check `bw_fct` is numeric and positive
+    if not isinstance(bw_fct, (int, float, np.integer, np.floating)):
+        raise TypeError(f"`bw_fct` must be a positive number, not an object of {type(bw_fct)}.")
 
-    # Assert `bw_fct` is numeric and positive
-    assert isinstance(bw_fct, (int, float))
-    assert bw_fct > 0
+    if not bw_fct > 0:
+        raise ValueError(f"`bw_fct` must be a positive number, not {bw_fct}.")
 
     # Preliminary calculations
     x_len = len(x)
@@ -675,22 +630,20 @@ def kde_linear(
     x_range = x_max - x_min
 
     # Determine grid
-    grid_min, grid_max, grid_len = get_grid(
+    grid_min, grid_max, grid_len = _get_grid(
         x_min, x_max, x_std, extend_fct, grid_len,
         custom_lims, extend, bound_correction
     )
     grid_counts, _, grid_edges = histogram(x, grid_len, (grid_min, grid_max))
 
     # Bandwidth estimation
-    bw = bw_fct * get_bw(x, bw, grid_counts, x_std, x_range)
+    bw = bw_fct * _get_bw(x, bw, grid_counts, x_std, x_range)
 
     # Density estimation
     if adaptive:
-        grid, pdf = _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len,
-                                  bound_correction)
+        grid, pdf = _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len,bound_correction)
     else:
-        grid, pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len,
-                                     bound_correction)
+        grid, pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
 
     if cumulative:
         pdf = pdf.cumsum() / pdf.sum()
@@ -701,14 +654,15 @@ def kde_linear(
         return grid, pdf
 
 
-def kde_circular(
+def _kde_circular(
     x,
     bw="taylor",
     bw_fct=1,
     bw_return=False,
     custom_lims=None,
     cumulative=False,
-    grid_len=512
+    grid_len=512,
+    **kwargs
 ):
     """
     1 dimensional density estimation for circular data.
@@ -746,14 +700,17 @@ def kde_circular(
         Defaults to 512.
     """
 
-    x = check_type(x)
+    x = _check_type(x)
 
     # All values between -pi and pi
-    x = normalize_angle(x)
+    x = numeric_utils._normalize_angle(x)
+
+    # Check `bw_fct` is numeric and positive
+    if not isinstance(bw_fct, (int, float, np.integer, np.floating)):
+        raise TypeError(f"`bw_fct` must be a positive number, not an object of {type(bw_fct)}.")
 
-    # Assert `bw_fct` is numeric and positive
-    assert isinstance(bw_fct, (int, float))
-    assert bw_fct > 0
+    if not bw_fct > 0:
+        raise ValueError(f"`bw_fct` must be a positive number, not {bw_fct}.")
 
     # Determine bandwidth
     if isinstance(bw, bool):
@@ -766,7 +723,7 @@ def kde_circular(
             raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
     if isinstance(bw, str):
         if bw == "taylor":
-            bw = bw_taylor(x)
+            bw = _bw_taylor(x)
         else:
             raise ValueError((
                 f"`bw` must be a positive numeric or `taylor`, not {bw}"
@@ -775,7 +732,7 @@ def kde_circular(
 
     # Determine grid
     if custom_lims is not None:
-        custom_lims = check_custom_lims(custom_lims, x.min(), x.max())
+        custom_lims = _check_custom_lims(custom_lims, x.min(), x.max())
         grid_min = custom_lims[0]
         grid_max = custom_lims[1]
         assert grid_min >= -np.pi, "Lower limit can't be smaller than -pi"
@@ -788,7 +745,7 @@ def kde_circular(
     bin_counts, _, bin_edges = histogram(x, bins=bins)
     grid = 0.5 * (bin_edges[1:] + bin_edges[:-1])
 
-    kern = vonmises_pdf(x=grid, mu=0, kappa=bw)
+    kern = _vonmises_pdf(x=grid, mu=0, kappa=bw)
     pdf = np.fft.fftshift(np.fft.irfft(np.fft.rfft(kern) * np.fft.rfft(bin_counts)))
     pdf /= len(x)
 
@@ -805,7 +762,7 @@ def _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
     """
     1 dimensional Gaussian kernel density estimation via
     convolution of the binned relative frequencies and a Gaussian filter.
-    This is an internal function used by `kde()`.
+    This is an internal function used by `_kde()`.
     """
 
     # Calculate relative frequencies per bin
@@ -842,11 +799,11 @@ def _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
     The implementation uses the binning technique.
     Since there is not an unique `bw`, the convolution is not possible.
     The alternative implemented in this function is known as Abramson's method.
-    This is an internal function used by `kde()`.
+    This is an internal function used by `_kde()`.
     """
     # Pilot computations used for bandwidth adjustment
-    pilot_grid, pilot_pdf = _kde_convolution(x, bw, grid_edges, grid_counts,
-                                             grid_len, bound_correction)
+    pilot_grid, pilot_pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len,
+                                             bound_correction)
 
     # Adds to avoid np.log(0) and zero division
     pilot_pdf += 1e-9
diff --git a/arviz/numeric_utils.py b/arviz/numeric_utils.py
index 35f66266e5..8e8fc9f6f8 100644
--- a/arviz/numeric_utils.py
+++ b/arviz/numeric_utils.py
@@ -1,11 +1,24 @@
 """Numerical utility functions for ArviZ."""
 import warnings
 import numpy as np
-from scipy.signal import convolve, convolve2d
+from scipy.signal import convolve2d
 from scipy.sparse import coo_matrix
 
 from .stats.stats_utils import histogram
 from .utils import _stack, _dot, _cov
+from .kde_utils import _kde
+
+def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):
+    """Deprecated Kernel Density Estimate
+
+    """
+    grid, pdf = _kde(x)
+
+    warnings.warn(
+        "_fast_kde() has been replaced by _kde() in kde_utils.py",
+        DeprecationWarning
+    )
+    return grid, pdf
 
 def _fast_kde_2d(x, y, gridsize=(128, 128), circular=False):
     """
@@ -141,3 +154,29 @@ def _sturges_formula(dataset, mult=1):
         Number of bins to use
     """
     return int(np.ceil(mult * np.log2(dataset.draw.size)) + 1)
+
+
+def _circular_mean(x, na_rm=False):
+    """
+    Computes mean of circular variable measured in radians.
+    The result is between -pi and pi.
+    """
+    if na_rm:
+        x = x[~np.isnan(x)]
+
+    sinr = np.sum(np.sin(x))
+    cosr = np.sum(np.cos(x))
+    mean = np.arctan2(sinr, cosr)
+
+    return mean
+
+def _normalize_angle(x, zero_centered=True):
+    """
+    Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
+    depending on `zero_centered`.
+    """
+    if zero_centered:
+        return (x + np.pi) % (2 * np.pi) - np.pi
+    else:
+        return x % (2 * np.pi)
+        
diff --git a/arviz/plots/backends/bokeh/bpvplot.py b/arviz/plots/backends/bokeh/bpvplot.py
index d40560fd95..62d431fbbc 100644
--- a/arviz/plots/backends/bokeh/bpvplot.py
+++ b/arviz/plots/backends/bokeh/bpvplot.py
@@ -14,7 +14,7 @@
     is_valid_quantile,
     vectorized_to_hex,
 )
-from ....kde_utils import kde
+from ....kde_utils import _kde
 
 def plot_bpv(
     ax,
@@ -90,7 +90,7 @@ def plot_bpv(
 
         if kind == "p_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=-1)
-            x_s, tstat_pit_dens = kde(tstat_pit)
+            x_s, tstat_pit_dens = _kde(tstat_pit)
             ax_i.line(x_s, tstat_pit_dens, line_width=linewidth, line_color=color)
             # ax_i.set_yticks([])
             if reference is not None:
@@ -110,7 +110,7 @@ def plot_bpv(
 
         elif kind == "u_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=0)
-            x_s, tstat_pit_dens = kde(tstat_pit)
+            x_s, tstat_pit_dens = _kde(tstat_pit)
             ax_i.line(x_s, tstat_pit_dens, line_color=color)
             if reference is not None:
                 if reference == "analytical":
diff --git a/arviz/plots/backends/bokeh/densityplot.py b/arviz/plots/backends/bokeh/densityplot.py
index 6a50881e5a..75f6636b2f 100644
--- a/arviz/plots/backends/bokeh/densityplot.py
+++ b/arviz/plots/backends/bokeh/densityplot.py
@@ -17,7 +17,7 @@
 )
 from ....stats import hdi
 from ....numeric_utils import histogram, get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 
 def plot_density(
     ax,
@@ -25,6 +25,7 @@ def plot_density(
     to_plot,
     colors,
     bw,
+    circular,
     figsize,
     length_plotters,
     rows,
@@ -97,6 +98,7 @@ def plot_density(
                 label,
                 colors[m_idx],
                 bw,
+                circular,
                 line_width,
                 markersize,
                 hdi_prob,
@@ -124,6 +126,7 @@ def _d_helper(
     vname,
     color,
     bw,
+    circular,
     line_width,
     markersize,
     hdi_prob,
@@ -144,7 +147,7 @@ def _d_helper(
         else:
             new_vec = vec
 
-        x, density = kde(new_vec, bw_fct=bw)
+        x, density = _kde(new_vec, circular=circular, bw=bw)
         density *= hdi_prob
         xmin, xmax = x[0], x[-1]
         ymin, ymax = density[0], density[-1]
@@ -228,7 +231,7 @@ def _d_helper(
         plotted.append(ax.diamond(xmax, 0, line_color="black", fill_color=color, size=markersize))
 
     if point_estimate is not None:
-        est = calculate_point_estimate(point_estimate, vec, bw)
+        est = calculate_point_estimate(point_estimate, vec, bw, circular)
         plotted.append(ax.circle(est, 0, fill_color=color, line_color="black", size=markersize))
 
     _title = Title()
diff --git a/arviz/plots/backends/bokeh/distplot.py b/arviz/plots/backends/bokeh/distplot.py
index 8e5b0520da..bc851146e8 100644
--- a/arviz/plots/backends/bokeh/distplot.py
+++ b/arviz/plots/backends/bokeh/distplot.py
@@ -20,6 +20,7 @@ def plot_dist(
     rotated,
     rug,
     bw,
+    circular,
     quantiles,
     contour,
     fill_last,
@@ -89,6 +90,7 @@ def plot_dist(
             rug=rug,
             label=label,
             bw=bw,
+            circular=circular,
             quantiles=quantiles,
             rotated=rotated,
             contour=contour,
diff --git a/arviz/plots/backends/bokeh/forestplot.py b/arviz/plots/backends/bokeh/forestplot.py
index 9e151011d5..7649c9fcce 100644
--- a/arviz/plots/backends/bokeh/forestplot.py
+++ b/arviz/plots/backends/bokeh/forestplot.py
@@ -17,7 +17,7 @@
 from ....stats import hdi
 from ....stats.diagnostics import _ess, _rhat
 from ....numeric_utils import histogram, get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 from ....utils import conditional_jit
 
 
@@ -572,7 +572,7 @@ def ridgeplot(self, mult, ridgeplot_kind):
                 _, density, x = histogram(values, bins=bins)
                 x = x[:-1]
             elif kind == "density":
-                x, density = kde(values)
+                x, density = _kde(values)
                 density_q = density.cumsum() / density.sum()
 
             xvals.append(x)
diff --git a/arviz/plots/backends/bokeh/loopitplot.py b/arviz/plots/backends/bokeh/loopitplot.py
index 5eeec57d75..671afe3d3f 100644
--- a/arviz/plots/backends/bokeh/loopitplot.py
+++ b/arviz/plots/backends/bokeh/loopitplot.py
@@ -7,7 +7,7 @@
 
 from . import backend_kwarg_defaults
 from .. import show_layout
-from ....kde_utils import kde
+from ....kde_utils import _kde
 from ...plot_utils import _scale_fig_size
 
 def plot_loo_pit(
@@ -183,7 +183,7 @@ def plot_loo_pit(
             )
         else:
             for idx in range(n_unif):
-                x_s, unif_density = kde(unif[idx, :])
+                x_s, unif_density = _kde(unif[idx, :])
                 ax.line(
                     x_s,
                     unif_density,
diff --git a/arviz/plots/backends/bokeh/posteriorplot.py b/arviz/plots/backends/bokeh/posteriorplot.py
index b7f10960c9..e7a17f1c8e 100644
--- a/arviz/plots/backends/bokeh/posteriorplot.py
+++ b/arviz/plots/backends/bokeh/posteriorplot.py
@@ -28,6 +28,7 @@ def plot_posterior(
     figsize,
     plotters,
     bw,
+    circular,
     bins,
     kind,
     point_estimate,
@@ -75,6 +76,7 @@ def plot_posterior(
             selection,
             ax=ax_,
             bw=bw,
+            circular=circular,
             bins=bins,
             kind=kind,
             point_estimate=point_estimate,
@@ -104,6 +106,7 @@ def _plot_posterior_op(
     selection,
     ax,
     bw,
+    circular,
     linewidth,
     bins,
     kind,
@@ -192,7 +195,7 @@ def display_rope(max_data):
     def display_point_estimate(max_data):
         if not point_estimate:
             return
-        point_value = calculate_point_estimate(point_estimate, values, bw)
+        point_value = calculate_point_estimate(point_estimate, values, bw, circular)
         sig_figs = format_sig_figs(point_value, round_to)
         point_text = "{point_estimate}={point_value:.{sig_figs}g}".format(
             point_estimate=point_estimate, point_value=point_value, sig_figs=sig_figs
@@ -234,6 +237,7 @@ def format_axes():
         plot_kde(
             values,
             bw=bw,
+            circular=circular,
             fill_kwargs={"fill_alpha": kwargs.pop("fill_alpha", 0)},
             plot_kwargs=kwargs,
             ax=ax,
diff --git a/arviz/plots/backends/bokeh/ppcplot.py b/arviz/plots/backends/bokeh/ppcplot.py
index 34b58b31ed..d8c41cdc10 100644
--- a/arviz/plots/backends/bokeh/ppcplot.py
+++ b/arviz/plots/backends/bokeh/ppcplot.py
@@ -6,7 +6,7 @@
 from ...kdeplot import plot_kde
 from ...plot_utils import _create_axes_grid, _scale_fig_size
 from ....numeric_utils import histogram, get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 
 
 def plot_ppc(
@@ -94,7 +94,7 @@ def plot_ppc(
             for vals in pp_sampled_vals:
                 vals = np.array([vals]).flatten()
                 if dtype == "f":
-                    pp_x, pp_density = kde(vals)
+                    pp_x, pp_density = _kde(vals)
                     pp_densities.append(pp_density)
                     pp_xs.append(pp_x)
                 else:
diff --git a/arviz/plots/backends/bokeh/violinplot.py b/arviz/plots/backends/bokeh/violinplot.py
index fc58c1083f..e9004836e7 100644
--- a/arviz/plots/backends/bokeh/violinplot.py
+++ b/arviz/plots/backends/bokeh/violinplot.py
@@ -4,8 +4,8 @@
 
 from . import backend_kwarg_defaults
 from .. import show_layout
-from ....kde_utils import kde
-from ....numeric_utils import _fast_kde, histogram, get_bins
+from ....numeric_utils import histogram, get_bins
+from ....kde_utils import _kde
 from ...plot_utils import make_label, _create_axes_grid, _scale_fig_size
 from ....stats import hdi
 
@@ -24,6 +24,7 @@ def plot_violin(
     rug_kwargs,
     bw,
     textsize,
+    circular,
     hdi_prob,
     quartiles,
     backend_kwargs,
@@ -65,7 +66,7 @@ def plot_violin(
         if val[0].dtype.kind == "i":
             dens = cat_hist(val, rug, shade, ax_, **shade_kwargs)
         else:
-            dens = _violinplot(val, rug, shade, bw, ax_, **shade_kwargs)
+            dens = _violinplot(val, rug, shade, bw, circular, ax_, **shade_kwargs)
 
         if rug:
             rug_x = -np.abs(np.random.normal(scale=max(dens) / 3.5, size=len(val)))
@@ -100,9 +101,14 @@ def plot_violin(
     return ax
 
 
-def _violinplot(val, rug, shade, bw, ax, **shade_kwargs):
+def _violinplot(val, rug, shade, bw, circular, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    x, density = kde(val, bw_fct=bw)
+    if bw is "default":
+        if circular:
+            bw = "taylor"
+        else:
+            bw = "experimental"
+    x, density = _kde(val, circular=circular, bw=bw)
 
     if not rug:
         x = np.concatenate([x, x[::-1]])
diff --git a/arviz/plots/backends/matplotlib/bpvplot.py b/arviz/plots/backends/matplotlib/bpvplot.py
index e2b4c44e94..1523e5c09d 100644
--- a/arviz/plots/backends/matplotlib/bpvplot.py
+++ b/arviz/plots/backends/matplotlib/bpvplot.py
@@ -13,7 +13,7 @@
     is_valid_quantile,
     matplotlib_kwarg_dealiaser,
 )
-from ....kde_utils import kde
+from ....kde_utils import _kde
 
 def plot_bpv(
     ax,
@@ -76,7 +76,7 @@ def plot_bpv(
 
         if kind == "p_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=-1)
-            x_s, tstat_pit_dens = kde(tstat_pit)
+            x_s, tstat_pit_dens = _kde(tstat_pit)
             ax_i.plot(x_s, tstat_pit_dens, linewidth=linewidth, color=color)
             ax_i.set_yticks([])
             if reference is not None:
@@ -95,7 +95,7 @@ def plot_bpv(
 
         elif kind == "u_value":
             tstat_pit = np.mean(pp_vals <= obs_vals, axis=0)
-            x_s, tstat_pit_dens = kde(tstat_pit)
+            x_s, tstat_pit_dens = _kde(tstat_pit)
             ax_i.plot(x_s, tstat_pit_dens, color=color)
             if reference is not None:
                 if reference == "analytical":
diff --git a/arviz/plots/backends/matplotlib/densityplot.py b/arviz/plots/backends/matplotlib/densityplot.py
index e303046978..343ff5d6d6 100644
--- a/arviz/plots/backends/matplotlib/densityplot.py
+++ b/arviz/plots/backends/matplotlib/densityplot.py
@@ -13,15 +13,16 @@
     calculate_point_estimate,
 )
 from ....numeric_utils import get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 
 
 def plot_density(
-    ax,
+    ax, 
     all_labels,
     to_plot,
     colors,
     bw,
+    circular,
     figsize,
     length_plotters,
     rows,
@@ -75,6 +76,7 @@ def plot_density(
                 label,
                 colors[m_idx],
                 bw,
+                circular,
                 titlesize,
                 xt_labelsize,
                 linewidth,
@@ -102,7 +104,8 @@ def _d_helper(
     vec,
     vname,
     color,
-    bw_fct,
+    bw,
+    circular,
     titlesize,
     xt_labelsize,
     linewidth,
@@ -124,10 +127,11 @@ def _d_helper(
         variable name
     color : str
         matplotlib color
-    bw_fct: float, optional
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
     titlesize : float
         font size for title
     xt_labelsize : float
@@ -153,7 +157,7 @@ def _d_helper(
         else:
             new_vec = vec
 
-        x, density = kde(new_vec, bw_fct=bw_fct)
+        x, density = _kde(new_vec, circular=circular, bw=bw)
         density *= hdi_prob
         xmin, xmax = x[0], x[-1]
         ymin, ymax = density[0], density[-1]
@@ -179,7 +183,7 @@ def _d_helper(
         ax.plot(xmax, 0, hdi_markers, color=color, markeredgecolor="k", markersize=markersize)
 
     if point_estimate is not None:
-        est = calculate_point_estimate(point_estimate, vec, bw_fct)
+        est = calculate_point_estimate(point_estimate, vec, bw)
         ax.plot(est, 0, "o", color=color, markeredgecolor="k", markersize=markersize)
 
     ax.set_yticks([])
diff --git a/arviz/plots/backends/matplotlib/distplot.py b/arviz/plots/backends/matplotlib/distplot.py
index c6fe73cae2..bbcc064235 100644
--- a/arviz/plots/backends/matplotlib/distplot.py
+++ b/arviz/plots/backends/matplotlib/distplot.py
@@ -18,6 +18,7 @@ def plot_dist(
     rotated,
     rug,
     bw,
+    circular,
     quantiles,
     contour,
     fill_last,
@@ -81,6 +82,7 @@ def plot_dist(
             rug=rug,
             label=label,
             bw=bw,
+            circular=circular,
             quantiles=quantiles,
             rotated=rotated,
             contour=contour,
diff --git a/arviz/plots/backends/matplotlib/forestplot.py b/arviz/plots/backends/matplotlib/forestplot.py
index a2afaed95b..93a134f5cd 100644
--- a/arviz/plots/backends/matplotlib/forestplot.py
+++ b/arviz/plots/backends/matplotlib/forestplot.py
@@ -10,7 +10,7 @@
 from ....stats import hdi
 from ....stats.diagnostics import _ess, _rhat
 from ....numeric_utils import histogram, get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 from ...plot_utils import _scale_fig_size, xarray_var_iter, make_label
 from ....utils import conditional_jit
 
@@ -519,7 +519,7 @@ def ridgeplot(self, mult, ridgeplot_kind):
                 density_q = density.cumsum() / density.sum()
                 x = x[:-1]
             elif kind == "density":
-                x, density = kde(values)
+                x, density = _kde(values)
                 density_q = density.cumsum() / density.sum()
 
             xvals.append(x)
diff --git a/arviz/plots/backends/matplotlib/loopitplot.py b/arviz/plots/backends/matplotlib/loopitplot.py
index 9e62ec34f2..a6f442a929 100644
--- a/arviz/plots/backends/matplotlib/loopitplot.py
+++ b/arviz/plots/backends/matplotlib/loopitplot.py
@@ -10,8 +10,7 @@
     _scale_fig_size,
     matplotlib_kwarg_dealiaser,
 )
-from ....kde_utils import kde
-from ....numeric_utils import _fast_kde
+from ....kde_utils import _kde
 
 
 def plot_loo_pit(
@@ -119,7 +118,7 @@ def plot_loo_pit(
             ax.axhspan(*hdi_odds, **hdi_kwargs)
         else:
             for idx in range(n_unif):
-                x_s, unif_density = kde(unif[idx, :])
+                x_s, unif_density = _kde(unif[idx, :])
                 x_ss[idx] = x_s
                 u_dens[idx] = unif_density
             ax.plot(x_ss.T, u_dens.T, **plot_unif_kwargs)
diff --git a/arviz/plots/backends/matplotlib/posteriorplot.py b/arviz/plots/backends/matplotlib/posteriorplot.py
index 94820a7ffc..09f4004105 100644
--- a/arviz/plots/backends/matplotlib/posteriorplot.py
+++ b/arviz/plots/backends/matplotlib/posteriorplot.py
@@ -26,6 +26,7 @@ def plot_posterior(
     figsize,
     plotters,
     bw,
+    circular,
     bins,
     kind,
     point_estimate,
@@ -67,6 +68,7 @@ def plot_posterior(
             selection,
             ax=ax_,
             bw=bw,
+            circular=circular,
             bins=bins,
             kind=kind,
             point_estimate=point_estimate,
@@ -95,6 +97,7 @@ def _plot_posterior_op(
     selection,
     ax,
     bw,
+    circular,
     linewidth,
     bins,
     kind,
@@ -204,7 +207,7 @@ def display_rope():
     def display_point_estimate():
         if not point_estimate:
             return
-        point_value = calculate_point_estimate(point_estimate, values, bw)
+        point_value = calculate_point_estimate(point_estimate, values, bw, circular)
         sig_figs = format_sig_figs(point_value, round_to)
         point_text = "{point_estimate}={point_value:.{sig_figs}g}".format(
             point_estimate=point_estimate, point_value=point_value, sig_figs=sig_figs
@@ -269,6 +272,7 @@ def format_axes():
         plot_kde(
             values,
             bw=bw,
+            circular=circular,
             fill_kwargs={"alpha": kwargs.pop("fill_alpha", 0)},
             plot_kwargs=kwargs,
             ax=ax,
diff --git a/arviz/plots/backends/matplotlib/ppcplot.py b/arviz/plots/backends/matplotlib/ppcplot.py
index abedc3d7d3..03a64930f8 100644
--- a/arviz/plots/backends/matplotlib/ppcplot.py
+++ b/arviz/plots/backends/matplotlib/ppcplot.py
@@ -13,7 +13,7 @@
     _scale_fig_size,
 )
 from ....numeric_utils import histogram, get_bins
-from ....kde_utils import kde
+from ....kde_utils import _kde
 _log = logging.getLogger(__name__)
 
 
@@ -149,7 +149,7 @@ def plot_ppc(
             for vals in pp_sampled_vals:
                 vals = np.array([vals]).flatten()
                 if dtype == "f":
-                    pp_x, pp_density = kde(vals)
+                    pp_x, pp_density = _kde(vals)
                     pp_densities.append(pp_density)
                     pp_xs.append(pp_x)
                 else:
@@ -367,7 +367,7 @@ def _set_animation(
     if kind == "kde":
         length = len(pp_sampled_vals)
         if dtype == "f":
-            x_vals, y_vals = kde(pp_sampled_vals[0])
+            x_vals, y_vals = _kde(pp_sampled_vals[0])
             max_max = max([max(kde(pp_sampled_vals[i])[1]) for i in range(length)])
             ax.set_ylim(0, max_max)
             (line,) = ax.plot(x_vals, y_vals, **plot_kwargs)
diff --git a/arviz/plots/backends/matplotlib/violinplot.py b/arviz/plots/backends/matplotlib/violinplot.py
index 758aad9af4..6d9536dd6f 100644
--- a/arviz/plots/backends/matplotlib/violinplot.py
+++ b/arviz/plots/backends/matplotlib/violinplot.py
@@ -4,7 +4,7 @@
 
 from . import backend_show
 from ....stats import hdi
-from ....kde_utils import kde
+from ....kde_utils import _kde
 from ....numeric_utils import histogram, get_bins
 from ...plot_utils import make_label, _create_axes_grid, matplotlib_kwarg_dealiaser, _scale_fig_size
 
@@ -23,6 +23,7 @@ def plot_violin(
     rug_kwargs,
     bw,
     textsize,
+    circular,
     hdi_prob,
     quartiles,
     backend_kwargs,
@@ -60,7 +61,7 @@ def plot_violin(
         if val[0].dtype.kind == "i":
             dens = cat_hist(val, rug, shade, ax_, **shade_kwargs)
         else:
-            dens = _violinplot(val, rug, shade, bw, ax_, **shade_kwargs)
+            dens = _violinplot(val, rug, shade, bw, circular, ax_, **shade_kwargs)
 
         if rug:
             rug_x = -np.abs(np.random.normal(scale=max(dens) / 3.5, size=len(val)))
@@ -85,9 +86,14 @@ def plot_violin(
     return ax
 
 
-def _violinplot(val, rug, shade, bw, ax, **shade_kwargs):
+def _violinplot(val, rug, shade, bw, circular, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    x, density = kde(val, bw_fct=bw)
+    if bw is "default":
+        if circular:
+            bw = "taylor"
+        else:
+            bw = "experimental"
+    x, density = _kde(val, circular=circular, bw=bw)
 
     if not rug:
         x = np.concatenate([x, x[::-1]])
diff --git a/arviz/plots/densityplot.py b/arviz/plots/densityplot.py
index 2955e08745..8e2c6a85ea 100644
--- a/arviz/plots/densityplot.py
+++ b/arviz/plots/densityplot.py
@@ -25,7 +25,8 @@ def plot_density(
     outline=True,
     hdi_markers="",
     shade=0.0,
-    bw=1,
+    bw="default",
+    circular=False,
     figsize=None,
     textsize=None,
     ax=None,
@@ -77,10 +78,16 @@ def plot_density(
     shade : Optional[float]
         Alpha blending value for the shaded area under the curve, between 0 (no shade) and 1
         (opaque). Defaults to 0.
-    bw : Optional[float]
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: Optional[float or str]
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular
+        and "taylor" when it is.
+    circular: Optional[bool]
+        If True, it interprets the values passed are from a circular variable measured in radians
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False.
     figsize : Optional[Tuple[int, int]]
         Figure size. If None it will be defined automatically.
     textsize: Optional[float]
@@ -212,12 +219,19 @@ def plot_density(
         length_plotters = max_plots
     rows, cols = default_grid(length_plotters, max_cols=3)
 
+    if bw is "default":
+        if circular:
+            bw = "taylor"
+        else:
+            bw = "experimental"
+
     plot_density_kwargs = dict(
         ax=ax,
         all_labels=all_labels,
         to_plot=to_plot,
         colors=colors,
         bw=bw,
+        circular=circular,
         figsize=figsize,
         length_plotters=length_plotters,
         rows=rows,
diff --git a/arviz/plots/distplot.py b/arviz/plots/distplot.py
index 707a4ffb2b..f35b46021d 100644
--- a/arviz/plots/distplot.py
+++ b/arviz/plots/distplot.py
@@ -15,7 +15,8 @@ def plot_dist(
     label=None,
     rotated=False,
     rug=False,
-    bw=1,
+    bw="default",
+    circular=False,
     quantiles=None,
     contour=True,
     fill_last=True,
@@ -58,10 +59,16 @@ def plot_dist(
         Whether to rotate the 1D KDE plot 90 degrees.
     rug : bool
         If True adds a rugplot. Defaults to False. Ignored for 2D KDE
-    bw : float
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: Optional[float or str]
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular 
+        and "taylor" when it is.
+    circular: Optional[bool]
+        If True, it interprets the values passed are from a circular variable measured in radians
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False.
     quantiles : list
         Quantiles in ascending order used to segment the KDE. Use [.25, .5, .75] for quartiles.
         Defaults to None.
@@ -171,6 +178,7 @@ def plot_dist(
         rotated=rotated,
         rug=rug,
         bw=bw,
+        circular=circular,
         quantiles=quantiles,
         contour=contour,
         fill_last=fill_last,
diff --git a/arviz/plots/energyplot.py b/arviz/plots/energyplot.py
index 16df90e836..91da33dc03 100644
--- a/arviz/plots/energyplot.py
+++ b/arviz/plots/energyplot.py
@@ -12,7 +12,7 @@ def plot_energy(
     legend=True,
     fill_alpha=(1, 0.75),
     fill_color=("C0", "C5"),
-    bw=1,
+    bw="experimental",
     textsize=None,
     fill_kwargs=None,
     plot_kwargs=None,
@@ -43,10 +43,10 @@ def plot_energy(
     fill_color : tuple of valid matplotlib color
         Color for Marginal energy distribution and Energy transition distribution.
         Defaults to ('C0', 'C5')
-    bw : float
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental". Defaults to "experimental"
         Only works if `kind='kde'`
     textsize: float
         Text size scaling factor for labels, titles and lines. If None it will be autoscaled based
diff --git a/arviz/plots/kdeplot.py b/arviz/plots/kdeplot.py
index 7bfd21870b..f35ce1f1f6 100644
--- a/arviz/plots/kdeplot.py
+++ b/arviz/plots/kdeplot.py
@@ -4,7 +4,7 @@
 
 from ..data import InferenceData
 from ..numeric_utils import _fast_kde_2d
-from ..kde_utils import kde
+from ..kde_utils import _kde
 from .plot_utils import get_plotting_function
 from ..rcparams import rcParams
 
@@ -15,7 +15,9 @@ def plot_kde(
     cumulative=False,
     rug=False,
     label=None,
-    bw=1,
+    bw="default",
+    adaptive=False,
+    circular=False,
     quantiles=None,
     rotated=False,
     contour=True,
@@ -51,10 +53,19 @@ def plot_kde(
         If True adds a rugplot. Defaults to False. Ignored for 2D KDE
     label : string
         Text to include as part of the legend
-    bw: float, optional
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular 
+        and "taylor" when it is.
+    adaptive: bool, optional. 
+        If True, an adaptative bandwidth is used. Only valid for 1D KDE.
+        Defaults to False.
+    circular: bool, optional.
+        If True, it interprets `values` is a circular variable measured in radians
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False.
     quantiles : list
         Quantiles in ascending order used to segment the KDE.
         Use [.25, .5, .75] for quartiles. Defaults to None.
@@ -129,6 +140,35 @@ def plot_kde(
 
         >>> az.plot_kde(mu_posterior, rug=True)
 
+    Plot KDE with adaptive bandwidth
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_kde(mu_posterior, adaptive=True)
+
+    Plot KDE with a different bandwidth estimator
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_kde(mu_posterior, bw="scott")
+
+    Plot KDE with a bandwidth specified manually
+
+    .. plot::
+        :context: close-figs
+
+        >>> az.plot_kde(mu_posterior, bw=0.4)
+
+    Plot KDE for a circular variable
+
+    .. plot::
+        :context: close-figs
+
+        >>> rvs = np.random.vonmises(mu=np.pi, kappa=2, size=500)
+        >>> az.plot_kde(rvs, circular=True)
+
 
     Plot a cumulative distribution
 
@@ -180,6 +220,10 @@ def plot_kde(
         :context: close-figs
 
         >>> az.plot_kde(mu_posterior, values2=tau_posterior, contour=False)
+        
+    See Also
+    --------
+    plot_kde : Compute and plot a kernel density estimate.
 
     """
     if isinstance(values, xr.Dataset):
@@ -194,7 +238,14 @@ def plot_kde(
         )
 
     if values2 is None:
-        grid, density = kde(values, bw_fct=bw, cumulative=cumulative)
+        
+        if bw is "default":
+            if circular:
+                bw = "taylor"
+            else:
+                bw = "experimental"
+
+        grid, density = _kde(values, circular, bw=bw, adaptive=adaptive, cumulative=cumulative)
         lower, upper = grid[0], grid[-1]
 
         if cumulative:
diff --git a/arviz/plots/loopitplot.py b/arviz/plots/loopitplot.py
index 9d4823a042..395c90cdf1 100644
--- a/arviz/plots/loopitplot.py
+++ b/arviz/plots/loopitplot.py
@@ -3,7 +3,7 @@
 import scipy.stats as stats
 
 from ..stats import loo_pit as _loo_pit
-from ..kde_utils import kde
+from ..kde_utils import _kde
 from .plot_utils import get_plotting_function
 
 from ..rcparams import rcParams
@@ -161,7 +161,7 @@ def plot_loo_pit(
         )
         unif_ecdf = unif_ecdf / n_data_points
     else:
-        x_vals, loo_pit_kde = kde(loo_pit)
+        x_vals, loo_pit_kde = _kde(loo_pit)
 
         unif = np.random.uniform(size=(n_unif, loo_pit.size))
         if use_hdi:
diff --git a/arviz/plots/plot_utils.py b/arviz/plots/plot_utils.py
index 6aad9984f8..6a0630204a 100644
--- a/arviz/plots/plot_utils.py
+++ b/arviz/plots/plot_utils.py
@@ -13,7 +13,7 @@
 from scipy.stats import mode
 import xarray as xr
 
-from ..kde_utils import kde
+from ..kde_utils import _kde
 from ..rcparams import rcParams
 
 KwargSpec = Dict[str, Any]
@@ -611,7 +611,7 @@ def get_plotting_function(plot_name, plot_module, backend):
     return plotting_method
 
 
-def calculate_point_estimate(point_estimate, values, bw=1):
+def calculate_point_estimate(point_estimate, values, bw="default", circular=False):
     """Validate and calculate the point estimate.
 
     Parameters
@@ -620,10 +620,16 @@ def calculate_point_estimate(point_estimate, values, bw=1):
         Plot point estimate per variable. Values should be 'mean', 'median', 'mode' or None.
         Defaults to 'auto' i.e. it falls back to default set in rcParams.
     values : 1-d array
-    bw: float, optional
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: Optional[float or str]
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular 
+        and "taylor" when it is.
+    circular: Optional[bool]
+        If True, it interprets the values passed are from a circular variable measured in radians
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False.
 
     Returns
     -------
@@ -643,7 +649,12 @@ def calculate_point_estimate(point_estimate, values, bw=1):
         point_value = values.mean()
     elif point_estimate == "mode":
         if isinstance(values[0], float):
-            x, density = kde(values, bw_fct=bw)
+            if bw is "default":
+                if circular:
+                    bw = "taylor"
+                else:
+                    bw = "experimental"
+            x, density = _kde(values, circular=circular, bw=bw)
             point_value = x[np.argmax(density)]
         else:
             point_value = mode(values)[0][0]
@@ -712,7 +723,7 @@ def sample_reference_distribution(dist, shape):
     densities = []
     dist_rvs = dist.rvs(size=shape)
     for idx in range(shape[1]):
-        x_s, density = kde(dist_rvs[:, idx])
+        x_s, density = _kde(dist_rvs[:, idx])
         x_ss.append(x_s)
         densities.append(density)
     return np.array(x_ss).T, np.array(densities).T
diff --git a/arviz/plots/posteriorplot.py b/arviz/plots/posteriorplot.py
index d502c07d34..a2ddd9073b 100644
--- a/arviz/plots/posteriorplot.py
+++ b/arviz/plots/posteriorplot.py
@@ -26,7 +26,8 @@ def plot_posterior(
     rope=None,
     ref_val=None,
     kind="kde",
-    bw=1,
+    bw="default",
+    circular=False,
     bins=None,
     ax=None,
     backend=None,
@@ -82,10 +83,16 @@ def plot_posterior(
     kind: str
         Type of plot to display (kde or hist) For discrete variables this argument is ignored and
         a histogram is always used.
-    bw : float
-        Bandwidth scaling factor for 1D KDE. Must be larger than 0.
-        The higher this number the smoother the KDE will be.
-        Defaults to 1 which means the bandwidth is not modified.
+    bw: float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular 
+        and "taylor" when it is. Only works if `kind == kde`.
+    circular: bool, optional
+        If True, it interprets the values passed are from a circular variable measured in radians
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False. 
         Only works if `kind == kde`.
     bins: integer or sequence or 'auto', optional
         Controls the number of bins, accepts the same keywords `matplotlib.hist()` does. Only works
@@ -225,6 +232,7 @@ def plot_posterior(
         figsize=figsize,
         plotters=plotters,
         bw=bw,
+        circular=circular,
         bins=bins,
         kind=kind,
         point_estimate=point_estimate,
diff --git a/arviz/plots/violinplot.py b/arviz/plots/violinplot.py
index 81fa9a868c..5037a3bec0 100644
--- a/arviz/plots/violinplot.py
+++ b/arviz/plots/violinplot.py
@@ -19,7 +19,8 @@ def plot_violin(
     rug=False,
     hdi_prob=None,
     shade=0.35,
-    bw=1,
+    bw="default",
+    circular=False,
     sharex=True,
     sharey=True,
     figsize=None,
@@ -63,9 +64,16 @@ def plot_violin(
     shade: float
         Alpha blending value for the shaded area under the curve, between 0
         (no shade) and 1 (opaque). Defaults to 0
-    bw: float
-        Bandwidth scaling factor. Must be larger than 0. The higher this number the smoother the
-        KDE will be. Defaults to 1 (no modification)
+    bw: float or str, optional
+        If numeric, indicates the bandwidth and must be positive.
+        If str, indicates the method to estimate the bandwidth and must be
+        one of "scott", "silverman", "isj" or "experimental" when `circular` is False
+        and "taylor" (for now) when `circular` is True.
+        Defaults to "default" which means "experimental" when variable is not circular
+        and "taylor" when it is.
+    circular: bool, optional.
+        If True, it interprets `values` is a circular variable measured in radians
+        and a circular KDE is used. Defaults to False.
     figsize: tuple
         Figure size. If None it will be defined automatically.
     textsize: int
@@ -150,6 +158,7 @@ def plot_violin(
         rug_kwargs=rug_kwargs,
         bw=bw,
         textsize=textsize,
+        circular=circular,
         hdi_prob=hdi_prob,
         quartiles=quartiles,
         backend_kwargs=backend_kwargs,
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index 94a7a1072d..c089b4bdc2 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -24,7 +24,7 @@
     get_log_likelihood as _get_log_likelihood,
 )
 from ..numeric_utils import histogram, get_bins
-from ..kde_utils import kde
+from ..kde_utils import _kde
 from ..utils import _var_names, Numba, _numba_var, get_coords, credible_interval_warning
 from ..rcparams import rcParams
 
@@ -546,7 +546,7 @@ def _hdi_multimodal(ary, hdi_prob, skipna, max_modes):
         ary = ary[~np.isnan(ary)]
 
     if ary.dtype.kind == "f":
-        bins, density = kde(ary)
+        bins, density = _kde(ary)
         lower, upper = bins[0], bins[-1]
         range_x = upper - lower
         dx = range_x / len(density)
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index b7209e4de7..73275ddec0 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -45,7 +45,7 @@
     plot_bpv,
 )
 from ...utils import _cov
-from ...kde_utils import kde
+from ...kde_utils import _kde
 
 rcParams["data.load"] = "eager"
 
@@ -929,7 +929,7 @@ def test_kde_scipy(limits):
     and the implementation in scipy
     """
     data = np.random.normal(0, 1, 10000)
-    grid, density_own = kde(data, custom_lims=limits)[1]
+    grid, density_own = _kde(data, custom_lims=limits)[1]
     density_sp = gaussian_kde(data).evaluate(grid)
     np.testing.assert_almost_equal(density_own.sum(), density_sp.sum(), 1)
 
@@ -940,7 +940,7 @@ def test_kde_cumulative(limits):
     Evaluates if last value of cumulative density is 1
     """
     data = np.random.normal(0, 1, 1000)
-    density = kde(data, custom_lims=limits, cumulative=True)[1]
+    density = _kde(data, custom_lims=limits, cumulative=True)[1]
     np.testing.assert_almost_equal(round(density[-1], 3), 1)
 
 @pytest.mark.parametrize(
diff --git a/arviz/tests/base_tests/test_utils.py b/arviz/tests/base_tests/test_utils.py
index 2d507d9121..149b960bd1 100644
--- a/arviz/tests/base_tests/test_utils.py
+++ b/arviz/tests/base_tests/test_utils.py
@@ -5,6 +5,7 @@
 from unittest.mock import Mock
 import numpy as np
 import pytest
+import scipy.stats as st
 
 from arviz.data.base import dict_to_dataset
 from ...utils import (
@@ -17,6 +18,7 @@
     _subset_list,
 )
 from ...data import load_arviz_data, from_dict
+from ...numeric_utils import _circular_mean, _normalize_angle
 
 
 @pytest.fixture(scope="session")
@@ -271,3 +273,23 @@ def test_flatten_inference_data_to_dict(
                 assert not any("sample_stats" in item for item in res_dict)
         else:
             assert not any("prior" in item for item in res_dict)
+
+
+@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, - 2 * np.pi, -10 * np.pi])
+def test_circular_mean_scipy(mean):
+    """Test our `_circular_mean()` function gives same result than Scipy version."""
+    rvs = st.vonmises.rvs(loc=mean, kappa=1, size=1000)
+    mean_az = _circular_mean(rvs)
+    mean_sp = st.circmean(rvs, low=-np.pi, high=np.pi)
+    np.testing.assert_almost_equal(mean_az, mean_sp)
+
+
+@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, - 2 * np.pi, -10 * np.pi])
+def test_normalize_angle(mean):
+    """Testing _normalize_angles() return values between expected bounds"""
+    rvs = st.vonmises.rvs(loc=mean, kappa=1, size=1000)
+    values = _normalize_angle(rvs, zero_centered=True)
+    assert ((-np.pi <= values) & (values <= np.pi)).all()
+
+    values = _normalize_angle(rvs, zero_centered=False)
+    assert ((0 <= values) & (values <= 2 * np.pi)).all()
\ No newline at end of file

From 819e9a615477beff180aca42f5a25922e187876d Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Wed, 29 Jul 2020 20:00:58 -0300
Subject: [PATCH 06/11] DeprecationWarning replaced by FutureWarning

---
 arviz/numeric_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/arviz/numeric_utils.py b/arviz/numeric_utils.py
index 8e8fc9f6f8..46c8f95e55 100644
--- a/arviz/numeric_utils.py
+++ b/arviz/numeric_utils.py
@@ -16,7 +16,7 @@ def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):
 
     warnings.warn(
         "_fast_kde() has been replaced by _kde() in kde_utils.py",
-        DeprecationWarning
+        FutureWarning
     )
     return grid, pdf
 

From 6d1539406382a713adb2858d9c19e26b7da4eaf8 Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Mon, 3 Aug 2020 15:29:26 -0300
Subject: [PATCH 07/11] Minor fixes * custom_lims in _kde() can be a tuple but
 then it tried to modify an inmutable object. Now it is converted to a list
 internally. * test_kde_scipy() had a slicing step that was not necessary an
 caused an error when unpacking * No error is raised when the data passed has
 lenght of 0 or 1. It raises a warning and returns np.nan

---
 arviz/kde_utils.py                            | 19 +++++++++++++++----
 arviz/plots/backends/matplotlib/ppcplot.py    |  4 ++--
 .../tests/base_tests/test_plots_matplotlib.py |  2 +-
 3 files changed, 18 insertions(+), 7 deletions(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 9ab599f64a..57ab22ddac 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -313,9 +313,11 @@ def _check_type(x):
         raise e
 
     x = x[np.isfinite(x)]
+
     if x.size == 0:
         raise ValueError("`x` does not contain any finite number.")
-
+    if x.size == 1:
+        raise ValueError("`x` is of length 1. Can't produce a KDE with only one data point.")
     return x
 
 
@@ -338,7 +340,7 @@ def _check_custom_lims(custom_lims, x_min, x_max):
             f"`custom_lims` must be a numeric list or tuple of length 2.\n"
             f"Not an object of {type(custom_lims)}."
         ))
-
+    
     if len(custom_lims) != 2:
         raise AttributeError(
             f"`len(custom_lims)` must be 2, not {len(custom_lims)}.")
@@ -348,6 +350,7 @@ def _check_custom_lims(custom_lims, x_min, x_max):
         raise TypeError(
             "Elements of `custom_lims` must be numeric or None, not bool.")
 
+    custom_lims = list(custom_lims) # convert to a mutable object
     if custom_lims[0] is None:
         custom_lims[0] = x_min
 
@@ -613,7 +616,11 @@ def _kde_linear(
     """
 
     # Check `x` is from appropiate type
-    x = _check_type(x)
+    try:
+        x = _check_type(x)
+    except ValueError as e:
+        warnings.warn('Something failed: ' + str(e))
+        return np.array([np.nan]), np.array([np.nan])
 
     # Check `bw_fct` is numeric and positive
     if not isinstance(bw_fct, (int, float, np.integer, np.floating)):
@@ -700,7 +707,11 @@ def _kde_circular(
         Defaults to 512.
     """
 
-    x = _check_type(x)
+    try:
+        x = _check_type(x)
+    except ValueError as e:
+        warnings.warn('Something failed: ' + str(e))
+        return np.array([np.nan]), np.array([np.nan])
 
     # All values between -pi and pi
     x = numeric_utils._normalize_angle(x)
diff --git a/arviz/plots/backends/matplotlib/ppcplot.py b/arviz/plots/backends/matplotlib/ppcplot.py
index 03a64930f8..84faf60b4b 100644
--- a/arviz/plots/backends/matplotlib/ppcplot.py
+++ b/arviz/plots/backends/matplotlib/ppcplot.py
@@ -368,12 +368,12 @@ def _set_animation(
         length = len(pp_sampled_vals)
         if dtype == "f":
             x_vals, y_vals = _kde(pp_sampled_vals[0])
-            max_max = max([max(kde(pp_sampled_vals[i])[1]) for i in range(length)])
+            max_max = max([max(_kde(pp_sampled_vals[i])[1]) for i in range(length)])
             ax.set_ylim(0, max_max)
             (line,) = ax.plot(x_vals, y_vals, **plot_kwargs)
 
             def animate(i):
-                x_vals, y_vals = kde(pp_sampled_vals[i])
+                x_vals, y_vals = _kde(pp_sampled_vals[i])
                 line.set_data(x_vals, y_vals)
                 return line
 
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index 73275ddec0..357ab14f99 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -929,7 +929,7 @@ def test_kde_scipy(limits):
     and the implementation in scipy
     """
     data = np.random.normal(0, 1, 10000)
-    grid, density_own = _kde(data, custom_lims=limits)[1]
+    grid, density_own = _kde(data, custom_lims=limits)
     density_sp = gaussian_kde(data).evaluate(grid)
     np.testing.assert_almost_equal(density_own.sum(), density_sp.sum(), 1)
 

From 4e5bb7605b7373d8e9bd5f7fd2514ccf15064787 Mon Sep 17 00:00:00 2001
From: Tomas Capretto <tomicapretto@gmail.com>
Date: Mon, 3 Aug 2020 20:02:51 -0300
Subject: [PATCH 08/11] Minor modifications regarding last PR comments

---
 arviz/kde_utils.py | 15 ++++-----------
 1 file changed, 4 insertions(+), 11 deletions(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 57ab22ddac..9802a89ae7 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -297,17 +297,12 @@ def _check_type(x):
         sample data from the variable for which a density estimate is desired.
 
     """
-    if not isinstance(x, (list, np.ndarray)):
-        raise ValueError((
-            f"Can't produce a density estimator for {type(x)}.\n"
-            f"Please input a list with numbers or numpy array."
-        ))
     # Will raise an error if `x` can't be casted to numeric or flattened to one dimension.
     try:
         x = np.asfarray(x).flatten()
     except Exception as e:
-        print(
-            "The following exception occurred while trying to convert `x`",
+        warnings.warn(
+            "The following exception occurred while trying to convert `x`" + 
             "to a 1 dimensional float array."
         )
         raise e
@@ -390,14 +385,14 @@ def _get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=Tru
         Number of bins
     custom_lims: tuple or list
         Custom limits for the domain of the density estimation.
-        Must be numeric of length 2.
+        Must be numeric of length 2. Overrides `extend`.
     extend: bool, optional
         Whether to extend the range of the data or not.
         Default is True.
     bound_correction: bool, optional
         Whether the density estimations performs boundary correction or not.
         This does not impacts directly in the output, but is used
-        to override `extend`.
+        to override `extend`. Overrides `extend`.
         Default is False.
 
     Returns
@@ -417,8 +412,6 @@ def _get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=Tru
     grid_len = int(grid_len)
 
     # Set up domain
-    # `custom_lims` overrides `extend`
-    # `bound_correction` overrides `extend`
     if custom_lims is not None:
         custom_lims = _check_custom_lims(custom_lims, x_min, x_max)
         grid_min = custom_lims[0]

From 01e61c28fa3f0b50c2ce6e47b4fd35a109830241 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1s=20Capretto?= <tomicapretto@gmail.com>
Date: Tue, 4 Aug 2020 15:47:02 -0300
Subject: [PATCH 09/11] Update arviz/kde_utils.py

Co-authored-by: Oriol Abril-Pla <oriol.abril.pla@gmail.com>
---
 arviz/kde_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 9802a89ae7..9ed842c378 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -302,7 +302,7 @@ def _check_type(x):
         x = np.asfarray(x).flatten()
     except Exception as e:
         warnings.warn(
-            "The following exception occurred while trying to convert `x`" + 
+            "The following exception occurred while trying to convert `x`"
             "to a 1 dimensional float array."
         )
         raise e

From 0ee8cf34716734770aa470a018dc9f93352443eb Mon Sep 17 00:00:00 2001
From: Ari Hartikainen <ahartikainen@users.noreply.github.com>
Date: Wed, 5 Aug 2020 01:32:35 +0300
Subject: [PATCH 10/11] Minor fixes

---
 arviz/kde_utils.py                            | 196 +++++++++---------
 arviz/numeric_utils.py                        |  22 +-
 arviz/plots/backends/bokeh/bpvplot.py         |   1 +
 arviz/plots/backends/bokeh/densityplot.py     |   1 +
 arviz/plots/backends/bokeh/forestplot.py      |   1 -
 arviz/plots/backends/bokeh/loopitplot.py      |   1 +
 arviz/plots/backends/bokeh/violinplot.py      |   2 +-
 arviz/plots/backends/matplotlib/bpvplot.py    |   1 +
 .../plots/backends/matplotlib/densityplot.py  |   2 +-
 arviz/plots/backends/matplotlib/ppcplot.py    |   1 +
 arviz/plots/backends/matplotlib/violinplot.py |   2 +-
 arviz/plots/densityplot.py                    |   2 +-
 arviz/plots/distplot.py                       |   2 +-
 arviz/plots/kdeplot.py                        |  10 +-
 arviz/plots/plot_utils.py                     |   4 +-
 arviz/plots/posteriorplot.py                  |   4 +-
 arviz/stats/stats.py                          |   3 +-
 .../tests/base_tests/test_plots_matplotlib.py |   1 +
 arviz/tests/base_tests/test_utils.py          |   6 +-
 19 files changed, 136 insertions(+), 126 deletions(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index 9ed842c378..c9787f7ecb 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -1,17 +1,17 @@
+# pylint: disable=invalid-name
 """Kernel density estimation functions for ArviZ."""
 
 import warnings
 import numpy as np
 from scipy.fftpack import fft
-from scipy.optimize import fsolve
 from scipy.optimize import brentq
 from scipy.signal import gaussian, convolve
-from scipy.special import ive
+from scipy.special import ive  # pylint: disable=no-name-in-module
 
 from .stats.stats_utils import histogram
-import arviz.numeric_utils as numeric_utils
 
-def _bw_scott(x, grid_counts=None, x_std=None, x_range=None):
+
+def _bw_scott(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
     """Scott's Rule.
 
     """
@@ -21,7 +21,7 @@ def _bw_scott(x, grid_counts=None, x_std=None, x_range=None):
     return bw
 
 
-def _bw_silverman(x, grid_counts=None, x_std=None, x_range=None):
+def _bw_silverman(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
     """Silverman's Rule.
 
     """
@@ -78,6 +78,7 @@ def _bw_isj(x, grid_counts=None, x_std=None, x_range=None):
     h = t ** 0.5 * x_range
     return h
 
+
 def _bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
     """Experimental bandwidth estimator.
 
@@ -113,11 +114,12 @@ def _bw_taylor(x):
 
 
 _BW_METHODS_LINEAR = {
-        "scott": _bw_scott,
-        "silverman": _bw_silverman,
-        "isj": _bw_isj,
-        "experimental": _bw_experimental
-    }
+    "scott": _bw_scott,
+    "silverman": _bw_silverman,
+    "isj": _bw_isj,
+    "experimental": _bw_experimental,
+}
+
 
 def _get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
     """
@@ -139,30 +141,36 @@ def _get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
         Bandwidth
     """
     if isinstance(bw, bool):
-        raise ValueError((
-            f"`bw` must not be of type `bool`.\n"
-            f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(_BW_METHODS_LINEAR.keys())}."))
+        raise ValueError(
+            (
+                "`bw` must not be of type `bool`.\n"
+                "Expected a positive numeric or one of the following strings:\n"
+                "{}."
+            ).format(list(_BW_METHODS_LINEAR.keys()))
+        )
     if isinstance(bw, (int, float)):
         if bw < 0:
-            raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
+            raise ValueError("Numeric `bw` must be positive.\nInput: {:.4f}.".format(bw))
     elif isinstance(bw, str):
         bw_lower = bw.lower()
 
         if bw_lower not in _BW_METHODS_LINEAR.keys():
-            raise ValueError((
-                f"Unrecognized bandwidth method.\n"
-                f"Input is: {method}.\n"
-                f"Expected one of: {list(_BW_METHODS_LINEAR.keys())}."
-            ))
+            raise ValueError(
+                (
+                    "Unrecognized bandwidth method.\n" "Input is: {}.\n" "Expected one of: {}."
+                ).format(bw_lower, list(_BW_METHODS_LINEAR.keys()))
+            )
 
         bw_fun = _BW_METHODS_LINEAR[bw_lower]
-        bw = bw_fun(x, grid_counts, x_std, x_range)           
+        bw = bw_fun(x, grid_counts, x_std, x_range)
     else:
-        raise ValueError((
-            f"Unrecognized `bw` argument.\n"
-            f"Expected a positive numeric or one of the following strings:\n"
-            f"{list(_BW_METHODS_LINEAR.keys())}."))
+        raise ValueError(
+            (
+                "Unrecognized `bw` argument.\n"
+                "Expected a positive numeric or one of the following strings:\n"
+                "{}."
+            ).format(list(_BW_METHODS_LINEAR.keys()))
+        )
     return bw
 
 
@@ -189,7 +197,9 @@ def _a1inv(x):
 
 
 def _kappa_mle(x):
-    mean = numeric_utils._circular_mean(x)
+    from numeric_utils import _circular_mean
+
+    mean = _circular_mean(x)
     kappa = _a1inv(np.mean(np.cos(x - mean)))
     return kappa
 
@@ -216,10 +226,7 @@ def _dct1d(x):
 
     x = np.concatenate((x[even_increasing], x[odd_decreasing]))
 
-    w_1k = np.r_[
-        1,
-        (2 * np.exp(-(0 + 1j) * (np.arange(1, x_len)) * np.pi / (2 * x_len)))
-    ]
+    w_1k = np.r_[1, (2 * np.exp(-(0 + 1j) * (np.arange(1, x_len)) * np.pi / (2 * x_len)))]
     output = np.real(w_1k * fft(x))
 
     return output
@@ -264,7 +271,7 @@ def _root(function, N, args, x):
     while not found:
         try:
             bw, res = brentq(function, 0, 0.01, args=args, full_output=True, disp=False)
-            found = True if res.converged else False
+            found = res.converged
         except ValueError:
             bw = 0
             tol *= 2.0
@@ -331,21 +338,21 @@ def _check_custom_lims(custom_lims, x_min, x_max):
 
     """
     if not isinstance(custom_lims, (list, tuple)):
-        raise TypeError((
-            f"`custom_lims` must be a numeric list or tuple of length 2.\n"
-            f"Not an object of {type(custom_lims)}."
-        ))
-    
+        raise TypeError(
+            (
+                "`custom_lims` must be a numeric list or tuple of length 2.\n"
+                "Not an object of {}."
+            ).format(type(custom_lims))
+        )
+
     if len(custom_lims) != 2:
-        raise AttributeError(
-            f"`len(custom_lims)` must be 2, not {len(custom_lims)}.")
+        raise AttributeError("`len(custom_lims)` must be 2, not {}.".format(len(custom_lims)))
 
     any_bool = any(isinstance(i, bool) for i in custom_lims)
     if any_bool:
-        raise TypeError(
-            "Elements of `custom_lims` must be numeric or None, not bool.")
+        raise TypeError("Elements of `custom_lims` must be numeric or None, not bool.")
 
-    custom_lims = list(custom_lims) # convert to a mutable object
+    custom_lims = list(custom_lims)  # convert to a mutable object
     if custom_lims[0] is None:
         custom_lims[0] = x_min
 
@@ -354,19 +361,19 @@ def _check_custom_lims(custom_lims, x_min, x_max):
 
     all_numeric = all(isinstance(i, (int, float, np.integer, np.float)) for i in custom_lims)
     if not all_numeric:
-        raise TypeError((
-            f"Elements of `custom_lims` must be numeric or None.\n"
-            f"At least one of them is not."
-        ))
+        raise TypeError(
+            ("Elements of `custom_lims` must be numeric or None.\n" "At least one of them is not.")
+        )
 
     if not custom_lims[0] < custom_lims[1]:
-        raise AttributeError(
-            f"`custom_lims[0]` must be smaller than `custom_lims[1]`.")
+        raise AttributeError("`custom_lims[0]` must be smaller than `custom_lims[1]`.")
 
     return custom_lims
 
-def _get_grid(x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True,
-              bound_correction=False):
+
+def _get_grid(
+    x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True, bound_correction=False
+):
     """
     Computes the grid that bins the data used to estimate the density function
 
@@ -550,7 +557,6 @@ def _kde_linear(
     custom_lims=None,
     cumulative=False,
     grid_len=512,
-    **kwargs
 ):
     """
     1 dimensional density estimation for linear data.
@@ -612,15 +618,17 @@ def _kde_linear(
     try:
         x = _check_type(x)
     except ValueError as e:
-        warnings.warn('Something failed: ' + str(e))
+        warnings.warn("Something failed: " + str(e))
         return np.array([np.nan]), np.array([np.nan])
 
     # Check `bw_fct` is numeric and positive
     if not isinstance(bw_fct, (int, float, np.integer, np.floating)):
-        raise TypeError(f"`bw_fct` must be a positive number, not an object of {type(bw_fct)}.")
+        raise TypeError(
+            "`bw_fct` must be a positive number, not an object of {}.".format(type(bw_fct))
+        )
 
-    if not bw_fct > 0:
-        raise ValueError(f"`bw_fct` must be a positive number, not {bw_fct}.")
+    if bw_fct <= 0:
+        raise ValueError("`bw_fct` must be a positive number, not {}.".format(bw_fct))
 
     # Preliminary calculations
     x_len = len(x)
@@ -631,8 +639,7 @@ def _kde_linear(
 
     # Determine grid
     grid_min, grid_max, grid_len = _get_grid(
-        x_min, x_max, x_std, extend_fct, grid_len,
-        custom_lims, extend, bound_correction
+        x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend, bound_correction
     )
     grid_counts, _, grid_edges = histogram(x, grid_len, (grid_min, grid_max))
 
@@ -641,7 +648,7 @@ def _kde_linear(
 
     # Density estimation
     if adaptive:
-        grid, pdf = _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len,bound_correction)
+        grid, pdf = _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
     else:
         grid, pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
 
@@ -655,14 +662,7 @@ def _kde_linear(
 
 
 def _kde_circular(
-    x,
-    bw="taylor",
-    bw_fct=1,
-    bw_return=False,
-    custom_lims=None,
-    cumulative=False,
-    grid_len=512,
-    **kwargs
+    x, bw="taylor", bw_fct=1, bw_return=False, custom_lims=None, cumulative=False, grid_len=512,
 ):
     """
     1 dimensional density estimation for circular data.
@@ -703,35 +703,36 @@ def _kde_circular(
     try:
         x = _check_type(x)
     except ValueError as e:
-        warnings.warn('Something failed: ' + str(e))
+        warnings.warn("Something failed: " + str(e))
         return np.array([np.nan]), np.array([np.nan])
 
     # All values between -pi and pi
-    x = numeric_utils._normalize_angle(x)
+    from numeric_utils import _normalize_angle
+
+    x = _normalize_angle(x)
 
     # Check `bw_fct` is numeric and positive
     if not isinstance(bw_fct, (int, float, np.integer, np.floating)):
-        raise TypeError(f"`bw_fct` must be a positive number, not an object of {type(bw_fct)}.")
+        raise TypeError(
+            "`bw_fct` must be a positive number, not an object of {}.".format(type(bw_fct))
+        )
 
-    if not bw_fct > 0:
-        raise ValueError(f"`bw_fct` must be a positive number, not {bw_fct}.")
+    if bw_fct <= 0:
+        raise ValueError("`bw_fct` must be a positive number, not {}.".format(bw_fct))
 
     # Determine bandwidth
     if isinstance(bw, bool):
-        raise ValueError((
-            "`bw` can't be of type `bool`.\n"
-            "Expected a positive numeric or 'taylor'"
-        ))
+        raise ValueError(
+            ("`bw` can't be of type `bool`.\n" "Expected a positive numeric or 'taylor'")
+        )
     if isinstance(bw, (int, float)):
         if bw < 0:
-            raise ValueError(f"Numeric `bw` must be positive.\nInput: {bw:.4f}.")
+            raise ValueError("Numeric `bw` must be positive.\nInput: {:.4f}.".format(bw))
     if isinstance(bw, str):
         if bw == "taylor":
             bw = _bw_taylor(x)
         else:
-            raise ValueError((
-                f"`bw` must be a positive numeric or `taylor`, not {bw}"
-            ))
+            raise ValueError(("`bw` must be a positive numeric or `taylor`, not {}".format(bw)))
     bw *= bw_fct
 
     # Determine grid
@@ -786,12 +787,12 @@ def _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
 
     if bound_correction:
         npad = int(grid_len / 5)
-        f = np.concatenate([f[npad - 1::-1], f, f[grid_len:grid_len - npad - 1:-1]])
-        pdf = convolve(f, kernel, mode="same", method="direct")[npad:npad + grid_len]
+        f = np.concatenate([f[npad - 1 :: -1], f, f[grid_len : grid_len - npad - 1 : -1]])
+        pdf = convolve(f, kernel, mode="same", method="direct")[npad : npad + grid_len]
         pdf /= bw * (2 * np.pi) ** 0.5
     else:
         pdf = convolve(f, kernel, mode="same", method="direct")
-        pdf /= (bw * (2 * np.pi) ** 0.5)
+        pdf /= bw * (2 * np.pi) ** 0.5
 
     grid = (grid_edges[1:] + grid_edges[:-1]) / 2
     return grid, pdf
@@ -806,8 +807,9 @@ def _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
     This is an internal function used by `_kde()`.
     """
     # Pilot computations used for bandwidth adjustment
-    pilot_grid, pilot_pdf = _kde_convolution(x, bw, grid_edges, grid_counts, grid_len,
-                                             bound_correction)
+    pilot_grid, pilot_pdf = _kde_convolution(
+        x, bw, grid_edges, grid_counts, grid_len, bound_correction
+    )
 
     # Adds to avoid np.log(0) and zero division
     pilot_pdf += 1e-9
@@ -830,27 +832,27 @@ def _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
         grid_padded = np.linspace(
             grid_edges[0] - grid_pad,
             grid_edges[grid_len - 1] + grid_pad,
-            num = grid_len + 2 * grid_npad
+            num=grid_len + 2 * grid_npad,
         )
-        grid_counts = np.concatenate([
-            grid_counts[grid_npad - 1:: -1],
-            grid_counts,
-            grid_counts[grid_len:grid_len - grid_npad - 1: -1]]
+        grid_counts = np.concatenate(
+            [
+                grid_counts[grid_npad - 1 :: -1],
+                grid_counts,
+                grid_counts[grid_len : grid_len - grid_npad - 1 : -1],
+            ]
         )
-        bw_adj = np.concatenate([
-            bw_adj[grid_npad - 1:: -1],
-            bw_adj,
-            bw_adj[grid_len:grid_len - grid_npad - 1: -1]]
+        bw_adj = np.concatenate(
+            [bw_adj[grid_npad - 1 :: -1], bw_adj, bw_adj[grid_len : grid_len - grid_npad - 1 : -1]]
         )
-        pdf_mat = ((grid_padded - grid_padded[:, None]) / bw_adj[:, None])
+        pdf_mat = (grid_padded - grid_padded[:, None]) / bw_adj[:, None]
         pdf_mat = np.exp(-0.5 * pdf_mat ** 2) * grid_counts[:, None]
-        pdf_mat /= ((2 * np.pi) ** 0.5 * bw_adj[:, None])
-        pdf = np.sum(pdf_mat[:, grid_npad:grid_npad + grid_len], axis=0) / len(x)
+        pdf_mat /= (2 * np.pi) ** 0.5 * bw_adj[:, None]
+        pdf = np.sum(pdf_mat[:, grid_npad : grid_npad + grid_len], axis=0) / len(x)
 
     else:
-        pdf_mat = ((grid - grid[:, None]) / bw_adj[:, None])
+        pdf_mat = (grid - grid[:, None]) / bw_adj[:, None]
         pdf_mat = np.exp(-0.5 * pdf_mat ** 2) * grid_counts[:, None]
-        pdf_mat /= ((2 * np.pi) ** 0.5 * bw_adj[:, None])
+        pdf_mat /= (2 * np.pi) ** 0.5 * bw_adj[:, None]
         pdf = np.sum(pdf_mat, axis=0) / len(x)
 
     return grid, pdf
diff --git a/arviz/numeric_utils.py b/arviz/numeric_utils.py
index 46c8f95e55..6eb88f180f 100644
--- a/arviz/numeric_utils.py
+++ b/arviz/numeric_utils.py
@@ -4,22 +4,25 @@
 from scipy.signal import convolve2d
 from scipy.sparse import coo_matrix
 
-from .stats.stats_utils import histogram
 from .utils import _stack, _dot, _cov
 from .kde_utils import _kde
+from .stats.stats_utils import histogram  # pylint: disable=unused-import
 
-def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):
+
+def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):  # pylint: disable=unused-argument
     """Deprecated Kernel Density Estimate
 
     """
-    grid, pdf = _kde(x)
+    if not (xmin is None and xmax is None):
+        custom_lims = (xmin, xmax)
+    else:
+        custom_lims = None
+    grid, pdf = _kde(x, cumulative=cumulative, bw=bw, custom_lims=custom_lims)
 
-    warnings.warn(
-        "_fast_kde() has been replaced by _kde() in kde_utils.py",
-        FutureWarning
-    )
+    warnings.warn("_fast_kde() has been replaced by _kde() in kde_utils.py", FutureWarning)
     return grid, pdf
 
+
 def _fast_kde_2d(x, y, gridsize=(128, 128), circular=False):
     """
     2D fft-based Gaussian kernel density estimate (KDE).
@@ -125,8 +128,7 @@ def get_bins(values):
     bins_sturges = (x_max - x_min) / (np.log2(values.size) + 1)
 
     # The Freedman-Diaconis histogram bin estimator.
-    iqr = np.subtract(
-        *np.percentile(values, [75, 25]))  # pylint: disable=assignment-from-no-return
+    iqr = np.subtract(*np.percentile(values, [75, 25]))  # pylint: disable=assignment-from-no-return
     bins_fd = 2 * iqr * values.size ** (-1 / 3)
 
     width = np.round(np.max([1, bins_sturges, bins_fd])).astype(int)
@@ -170,6 +172,7 @@ def _circular_mean(x, na_rm=False):
 
     return mean
 
+
 def _normalize_angle(x, zero_centered=True):
     """
     Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
@@ -179,4 +182,3 @@ def _normalize_angle(x, zero_centered=True):
         return (x + np.pi) % (2 * np.pi) - np.pi
     else:
         return x % (2 * np.pi)
-        
diff --git a/arviz/plots/backends/bokeh/bpvplot.py b/arviz/plots/backends/bokeh/bpvplot.py
index 62d431fbbc..75157de06e 100644
--- a/arviz/plots/backends/bokeh/bpvplot.py
+++ b/arviz/plots/backends/bokeh/bpvplot.py
@@ -16,6 +16,7 @@
 )
 from ....kde_utils import _kde
 
+
 def plot_bpv(
     ax,
     length_plotters,
diff --git a/arviz/plots/backends/bokeh/densityplot.py b/arviz/plots/backends/bokeh/densityplot.py
index 75f6636b2f..8910912dd1 100644
--- a/arviz/plots/backends/bokeh/densityplot.py
+++ b/arviz/plots/backends/bokeh/densityplot.py
@@ -19,6 +19,7 @@
 from ....numeric_utils import histogram, get_bins
 from ....kde_utils import _kde
 
+
 def plot_density(
     ax,
     all_labels,
diff --git a/arviz/plots/backends/bokeh/forestplot.py b/arviz/plots/backends/bokeh/forestplot.py
index 7649c9fcce..ce8c19323d 100644
--- a/arviz/plots/backends/bokeh/forestplot.py
+++ b/arviz/plots/backends/bokeh/forestplot.py
@@ -21,7 +21,6 @@
 from ....utils import conditional_jit
 
 
-
 def pairwise(iterable):
     """From itertools cookbook. [a, b, c, ...] -> (a, b), (b, c), ..."""
     first, second = tee(iterable)
diff --git a/arviz/plots/backends/bokeh/loopitplot.py b/arviz/plots/backends/bokeh/loopitplot.py
index 671afe3d3f..134e2765c4 100644
--- a/arviz/plots/backends/bokeh/loopitplot.py
+++ b/arviz/plots/backends/bokeh/loopitplot.py
@@ -10,6 +10,7 @@
 from ....kde_utils import _kde
 from ...plot_utils import _scale_fig_size
 
+
 def plot_loo_pit(
     ax,
     figsize,
diff --git a/arviz/plots/backends/bokeh/violinplot.py b/arviz/plots/backends/bokeh/violinplot.py
index e9004836e7..48ccfc3bb6 100644
--- a/arviz/plots/backends/bokeh/violinplot.py
+++ b/arviz/plots/backends/bokeh/violinplot.py
@@ -103,7 +103,7 @@ def plot_violin(
 
 def _violinplot(val, rug, shade, bw, circular, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    if bw is "default":
+    if bw == "default":
         if circular:
             bw = "taylor"
         else:
diff --git a/arviz/plots/backends/matplotlib/bpvplot.py b/arviz/plots/backends/matplotlib/bpvplot.py
index 1523e5c09d..08c695da7e 100644
--- a/arviz/plots/backends/matplotlib/bpvplot.py
+++ b/arviz/plots/backends/matplotlib/bpvplot.py
@@ -15,6 +15,7 @@
 )
 from ....kde_utils import _kde
 
+
 def plot_bpv(
     ax,
     length_plotters,
diff --git a/arviz/plots/backends/matplotlib/densityplot.py b/arviz/plots/backends/matplotlib/densityplot.py
index 343ff5d6d6..bfd803c248 100644
--- a/arviz/plots/backends/matplotlib/densityplot.py
+++ b/arviz/plots/backends/matplotlib/densityplot.py
@@ -17,7 +17,7 @@
 
 
 def plot_density(
-    ax, 
+    ax,
     all_labels,
     to_plot,
     colors,
diff --git a/arviz/plots/backends/matplotlib/ppcplot.py b/arviz/plots/backends/matplotlib/ppcplot.py
index 84faf60b4b..de5e060b38 100644
--- a/arviz/plots/backends/matplotlib/ppcplot.py
+++ b/arviz/plots/backends/matplotlib/ppcplot.py
@@ -14,6 +14,7 @@
 )
 from ....numeric_utils import histogram, get_bins
 from ....kde_utils import _kde
+
 _log = logging.getLogger(__name__)
 
 
diff --git a/arviz/plots/backends/matplotlib/violinplot.py b/arviz/plots/backends/matplotlib/violinplot.py
index 6d9536dd6f..e5e3a0b9ca 100644
--- a/arviz/plots/backends/matplotlib/violinplot.py
+++ b/arviz/plots/backends/matplotlib/violinplot.py
@@ -88,7 +88,7 @@ def plot_violin(
 
 def _violinplot(val, rug, shade, bw, circular, ax, **shade_kwargs):
     """Auxiliary function to plot violinplots."""
-    if bw is "default":
+    if bw == "default":
         if circular:
             bw = "taylor"
         else:
diff --git a/arviz/plots/densityplot.py b/arviz/plots/densityplot.py
index 8e2c6a85ea..3fdaba47c1 100644
--- a/arviz/plots/densityplot.py
+++ b/arviz/plots/densityplot.py
@@ -219,7 +219,7 @@ def plot_density(
         length_plotters = max_plots
     rows, cols = default_grid(length_plotters, max_cols=3)
 
-    if bw is "default":
+    if bw == "default":
         if circular:
             bw = "taylor"
         else:
diff --git a/arviz/plots/distplot.py b/arviz/plots/distplot.py
index f35b46021d..2472386cdc 100644
--- a/arviz/plots/distplot.py
+++ b/arviz/plots/distplot.py
@@ -64,7 +64,7 @@ def plot_dist(
         If str, indicates the method to estimate the bandwidth and must be
         one of "scott", "silverman", "isj" or "experimental" when `circular` is False
         and "taylor" (for now) when `circular` is True.
-        Defaults to "default" which means "experimental" when variable is not circular 
+        Defaults to "default" which means "experimental" when variable is not circular
         and "taylor" when it is.
     circular: Optional[bool]
         If True, it interprets the values passed are from a circular variable measured in radians
diff --git a/arviz/plots/kdeplot.py b/arviz/plots/kdeplot.py
index f35ce1f1f6..b65a2c86b1 100644
--- a/arviz/plots/kdeplot.py
+++ b/arviz/plots/kdeplot.py
@@ -58,9 +58,9 @@ def plot_kde(
         If str, indicates the method to estimate the bandwidth and must be
         one of "scott", "silverman", "isj" or "experimental" when `circular` is False
         and "taylor" (for now) when `circular` is True.
-        Defaults to "default" which means "experimental" when variable is not circular 
+        Defaults to "default" which means "experimental" when variable is not circular
         and "taylor" when it is.
-    adaptive: bool, optional. 
+    adaptive: bool, optional.
         If True, an adaptative bandwidth is used. Only valid for 1D KDE.
         Defaults to False.
     circular: bool, optional.
@@ -220,7 +220,7 @@ def plot_kde(
         :context: close-figs
 
         >>> az.plot_kde(mu_posterior, values2=tau_posterior, contour=False)
-        
+
     See Also
     --------
     plot_kde : Compute and plot a kernel density estimate.
@@ -238,8 +238,8 @@ def plot_kde(
         )
 
     if values2 is None:
-        
-        if bw is "default":
+
+        if bw == "default":
             if circular:
                 bw = "taylor"
             else:
diff --git a/arviz/plots/plot_utils.py b/arviz/plots/plot_utils.py
index 6a0630204a..f6818c12d8 100644
--- a/arviz/plots/plot_utils.py
+++ b/arviz/plots/plot_utils.py
@@ -625,7 +625,7 @@ def calculate_point_estimate(point_estimate, values, bw="default", circular=Fals
         If str, indicates the method to estimate the bandwidth and must be
         one of "scott", "silverman", "isj" or "experimental" when `circular` is False
         and "taylor" (for now) when `circular` is True.
-        Defaults to "default" which means "experimental" when variable is not circular 
+        Defaults to "default" which means "experimental" when variable is not circular
         and "taylor" when it is.
     circular: Optional[bool]
         If True, it interprets the values passed are from a circular variable measured in radians
@@ -649,7 +649,7 @@ def calculate_point_estimate(point_estimate, values, bw="default", circular=Fals
         point_value = values.mean()
     elif point_estimate == "mode":
         if isinstance(values[0], float):
-            if bw is "default":
+            if bw == "default":
                 if circular:
                     bw = "taylor"
                 else:
diff --git a/arviz/plots/posteriorplot.py b/arviz/plots/posteriorplot.py
index a2ddd9073b..6cf43c7d9f 100644
--- a/arviz/plots/posteriorplot.py
+++ b/arviz/plots/posteriorplot.py
@@ -88,11 +88,11 @@ def plot_posterior(
         If str, indicates the method to estimate the bandwidth and must be
         one of "scott", "silverman", "isj" or "experimental" when `circular` is False
         and "taylor" (for now) when `circular` is True.
-        Defaults to "default" which means "experimental" when variable is not circular 
+        Defaults to "default" which means "experimental" when variable is not circular
         and "taylor" when it is. Only works if `kind == kde`.
     circular: bool, optional
         If True, it interprets the values passed are from a circular variable measured in radians
-        and a circular KDE is used. Only valid for 1D KDE. Defaults to False. 
+        and a circular KDE is used. Only valid for 1D KDE. Defaults to False.
         Only works if `kind == kde`.
     bins: integer or sequence or 'auto', optional
         Controls the number of bins, accepts the same keywords `matplotlib.hist()` does. Only works
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index c089b4bdc2..c5b0370f4f 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -22,8 +22,9 @@
     stats_variance_2d as svar,
     _circular_standard_deviation,
     get_log_likelihood as _get_log_likelihood,
+    histogram,
 )
-from ..numeric_utils import histogram, get_bins
+from ..numeric_utils import get_bins
 from ..kde_utils import _kde
 from ..utils import _var_names, Numba, _numba_var, get_coords, credible_interval_warning
 from ..rcparams import rcParams
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index 357ab14f99..08f4e44139 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -943,6 +943,7 @@ def test_kde_cumulative(limits):
     density = _kde(data, custom_lims=limits, cumulative=True)[1]
     np.testing.assert_almost_equal(round(density[-1], 3), 1)
 
+
 @pytest.mark.parametrize(
     "kwargs",
     [
diff --git a/arviz/tests/base_tests/test_utils.py b/arviz/tests/base_tests/test_utils.py
index 149b960bd1..cef49000f6 100644
--- a/arviz/tests/base_tests/test_utils.py
+++ b/arviz/tests/base_tests/test_utils.py
@@ -275,7 +275,7 @@ def test_flatten_inference_data_to_dict(
             assert not any("prior" in item for item in res_dict)
 
 
-@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, - 2 * np.pi, -10 * np.pi])
+@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, -2 * np.pi, -10 * np.pi])
 def test_circular_mean_scipy(mean):
     """Test our `_circular_mean()` function gives same result than Scipy version."""
     rvs = st.vonmises.rvs(loc=mean, kappa=1, size=1000)
@@ -284,7 +284,7 @@ def test_circular_mean_scipy(mean):
     np.testing.assert_almost_equal(mean_az, mean_sp)
 
 
-@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, - 2 * np.pi, -10 * np.pi])
+@pytest.mark.parametrize("mean", [0, np.pi, 4 * np.pi, -2 * np.pi, -10 * np.pi])
 def test_normalize_angle(mean):
     """Testing _normalize_angles() return values between expected bounds"""
     rvs = st.vonmises.rvs(loc=mean, kappa=1, size=1000)
@@ -292,4 +292,4 @@ def test_normalize_angle(mean):
     assert ((-np.pi <= values) & (values <= np.pi)).all()
 
     values = _normalize_angle(rvs, zero_centered=False)
-    assert ((0 <= values) & (values <= 2 * np.pi)).all()
\ No newline at end of file
+    assert ((values >= 0) & (values <= 2 * np.pi)).all()

From c75989f505214080d1e139c1512cf90a38b01994 Mon Sep 17 00:00:00 2001
From: ahartikainen <ahartikainen@github.com>
Date: Thu, 6 Aug 2020 00:52:38 +0300
Subject: [PATCH 11/11] fix docstrings

---
 arviz/kde_utils.py     | 143 ++++++++++++++++++-----------------------
 arviz/numeric_utils.py |  15 ++---
 2 files changed, 70 insertions(+), 88 deletions(-)

diff --git a/arviz/kde_utils.py b/arviz/kde_utils.py
index c9787f7ecb..cc1428a9b8 100644
--- a/arviz/kde_utils.py
+++ b/arviz/kde_utils.py
@@ -12,9 +12,7 @@
 
 
 def _bw_scott(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
-    """Scott's Rule.
-
-    """
+    """Scott's Rule."""
     if x_std is None:
         x_std = np.std(x)
     bw = 1.06 * x_std * len(x) ** (-0.2)
@@ -22,9 +20,7 @@ def _bw_scott(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
 
 
 def _bw_silverman(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
-    """Silverman's Rule.
-
-    """
+    """Silverman's Rule."""
     if x_std is None:
         x_std = np.std(x)
     q75, q25 = np.percentile(x, [75, 25])
@@ -35,8 +31,7 @@ def _bw_silverman(x, x_std=None, **kwargs):  # pylint: disable=unused-argument
 
 
 def _bw_isj(x, grid_counts=None, x_std=None, x_range=None):
-    """
-    Improved Sheather-Jones bandwidth estimation.
+    """Improved Sheather-Jones bandwidth estimation.
 
     Improved Sheather and Jones method as explained in [1]_.
     This is an internal version pretended to be used by the KDE estimator.
@@ -49,7 +44,6 @@ def _bw_isj(x, grid_counts=None, x_std=None, x_range=None):
        Z. I. Botev, J. F. Grotowski, and D. P. Kroese.
        Ann. Statist. 38 (2010), no. 5, 2916--2957.
     """
-
     x_len = len(x)
     if x_range is None:
         x_min = np.min(x)
@@ -80,9 +74,7 @@ def _bw_isj(x, grid_counts=None, x_std=None, x_range=None):
 
 
 def _bw_experimental(x, grid_counts=None, x_std=None, x_range=None):
-    """Experimental bandwidth estimator.
-
-    """
+    """Experimental bandwidth estimator."""
     bw_silverman = _bw_silverman(x, x_std=x_std)
     bw_isj = _bw_isj(x, grid_counts=grid_counts, x_range=x_range)
     return 0.5 * (bw_silverman + bw_isj)
@@ -122,8 +114,8 @@ def _bw_taylor(x):
 
 
 def _get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
-    """
-    Computes bandwidth for a given data `x` and `bw`.
+    """Compute bandwidth for a given data `x` and `bw`.
+
     Also checks `bw` is correctly specified.
 
     Parameters
@@ -175,6 +167,7 @@ def _get_bw(x, bw, grid_counts=None, x_std=None, x_range=None):
 
 
 def _vonmises_pdf(x, mu, kappa):
+    """Calculate vonmises_pdf."""
     if kappa <= 0:
         raise ValueError("Argument 'kappa' must be positive.")
     pdf = 1 / (2 * np.pi * ive(0, kappa)) * np.exp(np.cos(x - mu) - 1) ** kappa
@@ -182,7 +175,8 @@ def _vonmises_pdf(x, mu, kappa):
 
 
 def _a1inv(x):
-    """
+    """Compute inverse function.
+
     Inverse function of the ratio of the first and
     zeroth order Bessel functions of the first kind.
 
@@ -205,8 +199,7 @@ def _kappa_mle(x):
 
 
 def _dct1d(x):
-    """
-    Discrete Cosine Transform in 1 Dimension
+    """Discrete Cosine Transform in 1 Dimension.
 
     Parameters
     ----------
@@ -218,7 +211,6 @@ def _dct1d(x):
     -------
     output : DTC transformed values
     """
-
     x_len = len(x)
 
     even_increasing = np.arange(0, x_len, 2)
@@ -233,8 +225,9 @@ def _dct1d(x):
 
 
 def _fixed_point(t, N, k_sq, a_sq):
-    """
-    Implementation of the function t-zeta*gamma^[l](t) derived from equation (30) in [1]
+    """Calculate t-zeta*gamma^[l](t).
+
+    Implementation of the function t-zeta*gamma^[l](t) derived from equation (30) in [1].
 
     References
     ----------
@@ -242,7 +235,6 @@ def _fixed_point(t, N, k_sq, a_sq):
        Z. I. Botev, J. F. Grotowski, and D. P. Kroese.
        Ann. Statist. 38 (2010), no. 5, 2916--2957.
     """
-
     k_sq = np.asfarray(k_sq, dtype=np.float64)
     a_sq = np.asfarray(a_sq, dtype=np.float64)
 
@@ -288,8 +280,8 @@ def _root(function, N, args, x):
 
 
 def _check_type(x):
-    """
-    Checks the input is of the correct type.
+    """Check the input is of the correct type.
+
     It only accepts numeric lists/numpy arrays of 1 dimension or something that
     can be flattened to 1 dimension.
 
@@ -302,7 +294,6 @@ def _check_type(x):
     x : 1-D numpy array
         If no error is thrown, a 1 dimensional array of
         sample data from the variable for which a density estimate is desired.
-
     """
     # Will raise an error if `x` can't be casted to numeric or flattened to one dimension.
     try:
@@ -324,8 +315,8 @@ def _check_type(x):
 
 
 def _check_custom_lims(custom_lims, x_min, x_max):
-    """
-    Checks whether `custom_lims` are of the correct type.
+    """Check if `custom_lims` are of the correct type.
+
     It accepts numeric lists/tuples of length 2.
 
     Parameters
@@ -335,7 +326,6 @@ def _check_custom_lims(custom_lims, x_min, x_max):
     Returns
     -------
     None: Object of type None
-
     """
     if not isinstance(custom_lims, (list, tuple)):
         raise TypeError(
@@ -374,8 +364,7 @@ def _check_custom_lims(custom_lims, x_min, x_max):
 def _get_grid(
     x_min, x_max, x_std, extend_fct, grid_len, custom_lims, extend=True, bound_correction=False
 ):
-    """
-    Computes the grid that bins the data used to estimate the density function
+    """Compute the grid that bins the data used to estimate the density function.
 
     Parameters
     ----------
@@ -410,9 +399,7 @@ def _get_grid(
         Minimum value of the grid
     grid_max: float
         Maximum value of the grid
-
     """
-
     # Set up number of bins.
     if grid_len < 100:
         grid_len = 100
@@ -434,7 +421,7 @@ def _get_grid(
 
 
 def _kde(x, circular=False, **kwargs):
-    """1 dimensional density estimation.
+    """One dimensional density estimation.
 
     It is a wrapper around `kde_linear()` and `kde_circular()`.
 
@@ -461,81 +448,80 @@ def _kde(x, circular=False, **kwargs):
     .. plot::
         :context: close-figs
 
-    	>>> import numpy as np
-	>>> import matplotlib.pyplot as plt
-	>>> from arviz.kde_utils import _kde
-	>>>
-	>>> rvs = np.random.gamma(shape=1.8, size=1000)
-	>>> grid, pdf = _kde(rvs)
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> import numpy as np
+    >>> import matplotlib.pyplot as plt
+    >>> from arviz.kde_utils import _kde
+    >>>
+    >>> rvs = np.random.gamma(shape=1.8, size=1000)
+    >>> grid, pdf = _kde(rvs)
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Density estimation for linear data with Silverman's rule bandwidth
     .. plot::
         :context: close-figs
 
-	>>> grid, pdf = _kde(rvs, bw="silverman")
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> grid, pdf = _kde(rvs, bw="silverman")
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Density estimation for linear data with scaled bandwidth
     .. plot::
         :context: close-figs
 
-	>>> # bw_fct > 1 means more smoothness.
-	>>> grid, pdf = _kde(rvs, bw_fct=2.5)
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> # bw_fct > 1 means more smoothness.
+    >>> grid, pdf = _kde(rvs, bw_fct=2.5)
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Default density estimation for linear data with extended limits
     .. plot::
         :context: close-figs
 
-	>>> grid, pdf = _kde(rvs, bound_correction=False, extend=True, extend_fct=0.5)
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> grid, pdf = _kde(rvs, bound_correction=False, extend=True, extend_fct=0.5)
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Default density estimation for linear data with custom limits
     .. plot::
         :context: close-figs
-	# It accepts tuples and lists of length 2.
-	>>> grid, pdf = __kde(rvs, bound_correction=False, custom_lims=(0, 10))
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    # It accepts tuples and lists of length 2.
+    >>> grid, pdf = __kde(rvs, bound_correction=False, custom_lims=(0, 10))
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Default density estimation for circular data
     .. plot::
         :context: close-figs
 
-	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
-	>>> grid, pdf = _kde(rvs, circular=True)
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
+    >>> grid, pdf = _kde(rvs, circular=True)
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Density estimation for circular data with scaled bandwidth
     .. plot::
         :context: close-figs
 
-	>>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
-	>>> # bw_fct > 1 means less smoothness.
-	>>> grid, pdf = _kde(rvs, circular=True, bw_fct=3)
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> rvs = np.random.vonmises(mu=np.pi, kappa=1, size=500)
+    >>> # bw_fct > 1 means less smoothness.
+    >>> grid, pdf = _kde(rvs, circular=True, bw_fct=3)
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     Density estimation for circular data with custom limits
     .. plot::
         :context: close-figs
-	>>> # This is still experimental, does not always work.
-	>>> rvs = np.random.vonmises(mu=0, kappa=30, size=500)
-	>>> grid, pdf = _kde(rvs, circular=True, custom_lims=(-1, 1))
-	>>> plt.plot(grid, pdf)
-	>>> plt.show()
+    >>> # This is still experimental, does not always work.
+    >>> rvs = np.random.vonmises(mu=0, kappa=30, size=500)
+    >>> grid, pdf = _kde(rvs, circular=True, custom_lims=(-1, 1))
+    >>> plt.plot(grid, pdf)
+    >>> plt.show()
 
     See Also
     --------
     plot_kde : Compute and plot a kernel density estimate.
     arviz.kde_utils._kde: Arviz KDE estimator
-
     """
     if circular:
         kde_fun = _kde_circular
@@ -558,8 +544,7 @@ def _kde_linear(
     cumulative=False,
     grid_len=512,
 ):
-    """
-    1 dimensional density estimation for linear data.
+    """One dimensional density estimation for linear data.
 
     Given an array of data points `x` it returns an estimate of
     the probability density function that generated the samples in `x`.
@@ -613,7 +598,6 @@ def _kde_linear(
     pdf : Numpy array for the density estimates.
     bw: optional, the estimated bandwidth.
     """
-
     # Check `x` is from appropiate type
     try:
         x = _check_type(x)
@@ -664,8 +648,7 @@ def _kde_linear(
 def _kde_circular(
     x, bw="taylor", bw_fct=1, bw_return=False, custom_lims=None, cumulative=False, grid_len=512,
 ):
-    """
-    1 dimensional density estimation for circular data.
+    """One dimensional density estimation for circular data.
 
     Given an array of data points `x` measured in radians,
     it returns an estimate of the probability density function that generated
@@ -699,7 +682,6 @@ def _kde_circular(
         (a.k.a. the length of the grid used in the estimation)
         Defaults to 512.
     """
-
     try:
         x = _check_type(x)
     except ValueError as e:
@@ -764,12 +746,12 @@ def _kde_circular(
 
 
 def _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
-    """
-    1 dimensional Gaussian kernel density estimation via
+    """Kernel density with convolution.
+
+    One dimensional Gaussian kernel density estimation via
     convolution of the binned relative frequencies and a Gaussian filter.
     This is an internal function used by `_kde()`.
     """
-
     # Calculate relative frequencies per bin
     bin_width = grid_edges[1] - grid_edges[0]
     f = grid_counts / bin_width / len(x)
@@ -799,8 +781,9 @@ def _kde_convolution(x, bw, grid_edges, grid_counts, grid_len, bound_correction)
 
 
 def _kde_adaptive(x, bw, grid_edges, grid_counts, grid_len, bound_correction):
-    """
-    1 dimensional adaptive Gaussian kernel density estimation.
+    """Compute Adaptive Kernel Density Estimation.
+
+    One dimensional adaptive Gaussian kernel density estimation.
     The implementation uses the binning technique.
     Since there is not an unique `bw`, the convolution is not possible.
     The alternative implemented in this function is known as Abramson's method.
diff --git a/arviz/numeric_utils.py b/arviz/numeric_utils.py
index 6eb88f180f..a007e614d4 100644
--- a/arviz/numeric_utils.py
+++ b/arviz/numeric_utils.py
@@ -10,9 +10,7 @@
 
 
 def _fast_kde(x, cumulative=False, bw=4.5, xmin=None, xmax=None):  # pylint: disable=unused-argument
-    """Deprecated Kernel Density Estimate
-
-    """
+    """Kernel Density Estimate, Deprecated."""
     if not (xmin is None and xmax is None):
         custom_lims = (xmin, xmax)
     else:
@@ -150,7 +148,7 @@ def _sturges_formula(dataset, mult=1):
     mult: float
         Used to scale the number of bins up or down. Default is 1 for Sturges' formula.
 
-     Returns
+    Returns
     -------
     int
         Number of bins to use
@@ -159,8 +157,8 @@ def _sturges_formula(dataset, mult=1):
 
 
 def _circular_mean(x, na_rm=False):
-    """
-    Computes mean of circular variable measured in radians.
+    """Compute mean of circular variable measured in radians.
+
     The result is between -pi and pi.
     """
     if na_rm:
@@ -174,8 +172,9 @@ def _circular_mean(x, na_rm=False):
 
 
 def _normalize_angle(x, zero_centered=True):
-    """
-    Takes angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
+    """Normalize angles.
+
+    Take angles in radians and normalize them to [-pi, pi) or [0, 2 * pi)
     depending on `zero_centered`.
     """
     if zero_centered: