def METHOD_NAME(x, method="KernSmooth", resolution=401):
    """**Bandwidth Selection for Density Estimation**
    Bandwidth selector for :func:`.density` estimation. See ``bw_method`` argument in
    :func:`.scipy.stats.gaussian_kde`.
    The ``"KernSmooth"`` method is adapted from the ``dpik()`` function from the *KernSmooth* R
    package. In this case, it estimates the optimal AMISE bandwidth using the direct plug-in method
    with 2 levels for the Parzen-Rosenblatt estimator with Gaussian kernel.
    Parameters
    -----------
    x : Union[list, np.array, pd.Series]
        A vector of values.
    method : float
        The bandwidth of the kernel. The larger the values, the smoother the estimation. Can be a
        number, or ``"scott"`` or ``"silverman"``
        (see ``bw_method`` argument in :func:`.scipy.stats.gaussian_kde`), or ``"KernSmooth"``.
    resolution : int
        Only when ``method="KernSmooth"``. The number of equally-spaced points over which binning
        is performed to obtain kernel functional approximation (see ``gridsize`` argument in ``KernSmooth::dpik()``).
    Returns
    -------
    float
        Bandwidth value.
    See Also
    --------
    density
    Examples
    --------
    .. ipython:: python
      import neurokit2 as nk
      x = np.random.normal(0, 1, size=100)
      bw = nk.density_bandwidth(x)
      bw
      nk.density_bandwidth(x, method="scott")
      nk.density_bandwidth(x, method=1)
      @savefig p_density_bandwidth.png scale=100%
      x, y = nk.density(signal, bandwidth=bw, show=True)
      @suppress
      plt.close()
    References
    ----------
    * Jones, W. M. (1995). Kernel Smoothing, Chapman & Hall.
    """
    if isinstance(method, str):
        method = method.lower()
    if isinstance(method, (float, int)) or method != "kernsmooth":
        return scipy.stats.gaussian_kde(x, bw_method=method).factor
    n = len(x)
    stdev = np.nanstd(x, ddof=1)
    iqr = np.diff(np.percentile(x, [25, 75]))[0] / 1.349
    scalest = min(stdev, iqr)
    data_scaled = (x - np.nanmean(x)) / scalest
    min_scaled = np.nanmin(data_scaled)
    max_scaled = np.nanmax(data_scaled)
    gcounts = _density_linearbinning(
        x=data_scaled,
        gpoints=np.linspace(min_scaled, max_scaled, resolution),
        truncate=True,
    )
    alpha = (2 * np.sqrt(2) ** 9 / (7 * n)) ** (1 / 9)
    psi6hat = _density_bkfe(gcounts, 6, alpha, min_scaled, max_scaled)
    alpha = (-3 * np.sqrt(2 / np.pi) / (psi6hat * n)) ** (1 / 7)
    psi4hat = _density_bkfe(gcounts, 4, alpha, min_scaled, max_scaled)
    delta_0 = 1 / ((4 * np.pi) ** (1 / 10))
    output = scalest * delta_0 * (1 / (psi4hat * n)) ** (1 / 5)
    return output