docs: def kde_multivariate(events_x, events_y, xout=None, yout=None, bw=None):
      """ Multivariate Kernel Density Estimation
  
      Parameters
      ----------
      events_x, events_y: 1D ndarray
          The input points for kernel density estimation. Input
          is flattened automatically.
      bw: tuple (bwx, bwy) or None
          The bandwith for kernel density estimation.
      xout, yout: ndarray
          The coordinates at which the KDE should be computed.
          If set to none, input coordinates are used.
  
      Returns
      -------
      density: ndarray, same shape as `xout`
          The KDE for the points in (xout, yout)
  
      See Also
      --------
      `statsmodels.nonparametric.kernel_density.KDEMultivariate`
      """
      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 bw is None:
          # divide by 2 to make it comparable to histogram KDE
          bw = (bin_width_doane(events_x) / 2,
                bin_width_doane(events_y) / 2)
  
      positions = np.vstack([xout.flatten(), yout.flatten()])
      estimator_ly = KDEMultivariate(data=[events_x.flatten(),
                                           events_y.flatten()],
                                     var_type='cc', bw=bw)
  
      density = estimator_ly.pdf(positions)
      return density.reshape(xout.shape)