title: How to combine two outputs of KernelDensity from sklearn (Python)
body: <p>I want to multiply two Kernel Density Estimates to get a composite likelihood ... The multiplication operator doesn't exist for sklearn.</p>
  <p>I have KDR for data from 3 independent sources and I want multiply KDE from each sources. The below code is from <a href="https://scikit-learn.org/stable/auto_examples/neighbors/plot_species_kde.html#sphx-glr-auto-examples-neighbors-plot-species-kde-py" rel="nofollow noreferrer">https://scikit-learn.org/stable/auto_examples/neighbors/plot_species_kde.html#sphx-glr-auto-examples-neighbors-plot-species-kde-py</a> and should run easily.</p>
  <p>Although I made some small modification to get get two KDEs and then multiply them which fails.</p>
  <p>Someone can help. Below code should run without error except the multiplication part.</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from sklearn.datasets import fetch_species_distributions
  from sklearn.neighbors import KernelDensity
  
  # if basemap is available, we'll use it.
  # otherwise, we'll improvise later...
  try:
      from mpl_toolkits.basemap import Basemap
      basemap = True
  except ImportError:
      basemap = False
  
  def construct_grids(batch):
      &quot;&quot;&quot;Construct the map grid from the batch object
  
      Parameters
      ----------
      batch : Batch object
          The object returned by :func:`fetch_species_distributions`
  
      Returns
      -------
      (xgrid, ygrid) : 1-D arrays
          The grid corresponding to the values in batch.coverages
      &quot;&quot;&quot;
      # x,y coordinates for corner cells
      xmin = batch.x_left_lower_corner + batch.grid_size
      xmax = xmin + (batch.Nx * batch.grid_size)
      ymin = batch.y_left_lower_corner + batch.grid_size
      ymax = ymin + (batch.Ny * batch.grid_size)
  
      # x coordinates of the grid cells
      xgrid = np.arange(xmin, xmax, batch.grid_size)
      # y coordinates of the grid cells
      ygrid = np.arange(ymin, ymax, batch.grid_size)
  
      return (xgrid, ygrid)
  
  # Get matrices/arrays of species IDs and locations
  data = fetch_species_distributions()
  species_names = ['Bradypus Variegatus', 'Microryzomys Minutus']
  
  Xtrain = np.vstack([data['train']['dd lat'], data['train']['dd long']]).T
  ytrain = np.array([d.decode('ascii').startswith('micro')
                    for d in data['train']['species']], dtype='int')
  Xtrain *= np.pi / 180.  # Convert lat/long to radians
  
  # Set up the data grid for the contour plot
  xgrid, ygrid = construct_grids(data)
  X, Y = np.meshgrid(xgrid[::5], ygrid[::5][::-1])
  land_reference = data.coverages[6][::5, ::5]
  land_mask = (land_reference &gt; -9999).ravel()
  
  xy = np.vstack([Y.ravel(), X.ravel()]).T
  xy = xy[land_mask]
  xy *= np.pi / 180.
  
  # Plot map of South America with distributions of each species
  fig = plt.figure()
  fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05)
  
  kde0 = KernelDensity(bandwidth=0.04, metric='haversine',
                      kernel='gaussian', algorithm='ball_tree')
  kde0.fit(Xtrain[ytrain == 0])
  
  kde1 = KernelDensity(bandwidth=0.04, metric='haversine',
                      kernel='gaussian', algorithm='ball_tree')
  kde1.fit(Xtrain[ytrain == 1])
  
  kde01=kde0*kde1
  
  plt.subplot(1, 1, 1)
  
  # evaluate only on the land: -9999 indicates ocean
  Z = np.full(land_mask.shape[0], -9999, dtype='int')
  Z[land_mask] = np.exp(kde01.score_samples(xy))
  Z = Z.reshape(X.shape)
  
  # plot contours of the density
  levels = np.linspace(0, Z.max(), 25)
  plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)
  
  plt.show()
  </code></pre>