docs: def kde_histogram(events_x, events_y, xout=None, yout=None, bins=None):
      """ Histogram-based Kernel Density Estimation
  
      Parameters
      ----------
      events_x, events_y: 1D ndarray
          The input points for kernel density estimation. Input
          is flattened automatically.
      xout, yout: ndarray
          The coordinates at which the KDE should be computed.
          If set to none, input coordinates are used.
      bins: tuple (binsx, binsy)
          The number of bins to use for the histogram.
  
      Returns
      -------
      density: ndarray, same shape as `xout`
          The KDE for the points in (xout, yout)
  
      See Also
      --------
      `numpy.histogram2d`
      `scipy.interpolate.RectBivariateSpline`
      """
      valid_combi = ((xout is None and yout is None) or
                     (xout is not None and yout is not None)
                     )
      if not valid_combi:
          raise ValueError("Both `xout` and `yout` must be (un)set.")
  
      if yout is None and yout is None:
          xout = events_x
          yout = events_y
  
      if bins is None:
          bins = (max(5, bin_num_doane(events_x)),
                  max(5, bin_num_doane(events_y)))
  
      # Compute the histogram
      hist2d, xedges, yedges = np.histogram2d(x=events_x,
                                              y=events_y,
                                              bins=bins,
                                              normed=True)
      xip = xedges[1:]-(xedges[1]-xedges[0])/2
      yip = yedges[1:]-(yedges[1]-yedges[0])/2
  
      estimator = RectBivariateSpline(x=xip, y=yip, z=hist2d)
      density = estimator.ev(xout, yout)
      density[density < 0] = 0
  
      return density.reshape(xout.shape)