docs: 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