title: Relation between 2D KDE bandwidth in sklearn vs bandwidth in scipy
body: <p>I'm attempting to compare the performance of <a href="http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html" rel="nofollow noreferrer">sklearn.neighbors.KernelDensity</a> versus <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html" rel="nofollow noreferrer">scipy.stats.gaussian_kde</a> for a two dimensional array.</p>
  
  <p>From <a href="http://jakevdp.github.io/blog/2013/12/01/kernel-density-estimation/" rel="nofollow noreferrer">this article</a> I see that the bandwidths (bw) are treated differently in each function. The article gives a recipe for setting the correct bw in <code>scipy</code> so it will be equivalent to the one used in <code>sklearn</code> . Basically it divides the bw by the sample standard deviation. The result is this:</p>
  
  <pre><code># For sklearn
  bw = 0.15
  
  # For scipy
  bw = 0.15/x.std(ddof=1)
  </code></pre>
  
  <p>where <code>x</code> is the sample array I'm using to obtain the KDE. This works just fine in 1D, but I can't make it work in 2D.</p>
  
  <p>Here's a <code>MWE</code> of what I got:</p>
  
  <pre><code>import numpy as np
  from scipy import stats
  from sklearn.neighbors import KernelDensity
  
  # Generate random data.
  n = 1000
  m1, m2 = np.random.normal(0.2, 0.2, size=n), np.random.normal(0.2, 0.2, size=n)
  # Define limits.
  xmin, xmax = min(m1), max(m1)
  ymin, ymax = min(m2), max(m2)
  # Format data.
  x, y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
  positions = np.vstack([x.ravel(), y.ravel()])
  values = np.vstack([m1, m2])
  
  # Define some point to evaluate the KDEs.
  x1, y1 = 0.5, 0.5
  
  # -------------------------------------------------------
  # Perform a kernel density estimate on the data using scipy.
  kernel = stats.gaussian_kde(values, bw_method=0.15/np.asarray(values).std(ddof=1))
  # Get KDE value for the point.
  iso1 = kernel((x1,y1))
  print 'iso1 = ', iso[0]
  
  # -------------------------------------------------------
  # Perform a kernel density estimate on the data using sklearn.
  kernel_sk = KernelDensity(kernel='gaussian', bandwidth=0.15).fit(zip(*values))
  # Get KDE value for the point.
  iso2 = kernel_sk.score_samples([[x1, y1]])
  print 'iso2 = ', np.exp(iso2[0])
  </code></pre>
  
  <p>( <code>iso2</code> is presented as an exponential since <code>sklearn</code> returns the log values)</p>
  
  <p>The results I get for <code>iso1</code> and <code>iso2</code> are different and I'm lost as to how should I affect the bandwidth (in either function) to make them equal (as they should).</p>
  
  <hr>
  
  <p><strong>Add</strong></p>
  
  <p>I was advised over at <code>sklearn</code> chat (by ep) that I should scale the values in <code>(x,y)</code> before calculating the kernel with <code>scipy</code> in order to obtain comparable results with <code>sklearn</code>.</p>
  
  <p>So this is what I did:</p>
  
  <pre><code># Scale values.
  x_val_sca = np.asarray(values[0])/np.asarray(values).std(axis=1)[0]
  y_val_sca = np.asarray(values[1])/np.asarray(values).std(axis=1)[1]
  values = [x_val_sca, y_val_sca]
  kernel = stats.gaussian_kde(values, bw_method=bw_value)
  </code></pre>
  
  <p>ie: I scaled both dimensions before getting the kernel with <code>scipy</code> while leaving the line that obtains the kernel in <code>sklearn</code> untouched.</p>
  
  <p>This gave better results but there's still differences in the kernels obtained:</p>
  
  <p><img src="https://i.stack.imgur.com/tmY3w.png" alt="kernels"></p>
  
  <p>where the red dot is the <code>(x1,y1)</code> point in the code. So as can be seen, there are still differences in the shapes of the density estimates, albeit very small ones. Perhaps this is the best that can be achieved?</p>