def _find_quantile_level(density, x, y, xp, yp, quantile, acc=.01,
                         ret_err=False):
    """Find density level for a given data quantile by iteration

    Parameters
    ----------
    density: 2d ndarray of shape (M, N)
        Kernel density estimate for which to compute the contours
    x: 2d ndarray of shape (M, N) or 1d ndarray of size M
        X-values corresponding to `kde`
    y: 2d ndarray of shape (M, N) or 1d ndarray of size M
        Y-values corresponding to `kde`
    xp: 1d ndarray of size D
        Event x-data from which to compute the quantile
    yp: 1d ndarray of size D
        Event y-data from which to compute the quantile
    quantile: float between 0 and 1
        Quantile along which to find contours in `kde` relative
        to its maximum
    acc: float
        Desired absolute accuracy (stopping criterion) of the
        contours
    ret_err: bool
        If True, also return the absolute error

    Returns
    -------
    level: float
        Contours level corresponding to the given quantile

    Notes
    -----
    A much more faster method (using interpolation) is implemented in
    :func:`get_quantile_levels`.
    NaN-values events in `xp` and `yp` are ignored.

    See Also
    --------
    skimage.measure.find_contours: Contour finding algorithm
    """
    if quantile >= 1 or quantile <= 0:
        raise ValueError("Invalid value for `quantile`: {}".format(quantile))

    # remove bad events
    bad = get_bad_vals(xp, yp)
    xp = xp[~bad]
    yp = yp[~bad]

    # initial guess
    level = quantile
    # error of current iteration
    err = 1
    # iteration factor (guarantees convergence)
    itfac = 1
    # total number of events
    nev = xp.size

    while np.abs(err) > acc:
        # compute contours
        conts = find_contours_level(density, x, y, level, closed=True)
        # compute number of points in contour
        isin = 0
        for ii in range(nev):
            for cc in conts:
                isin += PolygonFilter.point_in_poly((xp[ii], yp[ii]),
                                                    poly=cc)
                break  # no need to check other contours
        err = quantile - (nev - isin) / nev
        level += err * itfac
        itfac *= .9

    if ret_err:
        return level, err
    else:
        return level
Find density level for a given data quantile by iteration

    Parameters
    ----------
    density: 2d ndarray of shape (M, N)
        Kernel density estimate for which to compute the contours
    x: 2d ndarray of shape (M, N) or 1d ndarray of size M
        X-values corresponding to `kde`
    y: 2d ndarray of shape (M, N) or 1d ndarray of size M
        Y-values corresponding to `kde`
    xp: 1d ndarray of size D
        Event x-data from which to compute the quantile
    yp: 1d ndarray of size D
        Event y-data from which to compute the quantile
    quantile: float between 0 and 1
        Quantile along which to find contours in `kde` relative
        to its maximum
    acc: float
        Desired absolute accuracy (stopping criterion) of the
        contours
    ret_err: bool
        If True, also return the absolute error

    Returns
    -------
    level: float
        Contours level corresponding to the given quantile

    Notes
    -----
    A much more faster method (using interpolation) is implemented in
    :func:`get_quantile_levels`.
    NaN-values events in `xp` and `yp` are ignored.

    See Also
    --------
    skimage.measure.find_contours: Contour finding algorithm