SciPy KDE gradient

I am using the SciPy implementation of the kernel density estimate (KDE) (http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html), which is working fine so far. However, I would now like to obtain the gradient of the KDE at a particular point.

I have looked at the Python source for the library, but haven't been able to figure out how to easily implement this functionality. Is anybody aware of a method to do this?


If you look at the source you referenced, you'll see that the density estimation is constructed from contributions from all points in the dataset Assuming there's only one point points[:,i] you want to evaluate for the moment (lines 219–222):

diff = self.dataset - points[:, i, newaxis] tdiff = dot(self.inv_cov, diff) energy = sum(diff * tdiff, axis=0) / 2.0 result[i] = sum(exp(-energy), axis=0)

In matrix notation (no LaTeX available?), this would be written, for a single point D from the dataset and point p to be evaluated as

d = D - p t = Cov^-1 d e = 1/2 d^T t r = exp(-e)

The gradient you're looking for is grad(r) = (dr/dx, dr/dy):

dr/dx = d(exp(-e))/dx = -de/dx exp(-e) = -d(1/2 d^T Cov^-1 d)/dx exp(-e) = -(Cov^-1 d) exp(-e)

Likewise for dr/dy. Hence all you need to do is calculate the term Cov^-1 d and multiply it with the result you already obtained.

result = zeros((self.d,m), dtype=float) [...] diff = self.dataset - points[:, i, newaxis] tdiff = dot(self.inv_cov, diff) energy = sum(diff * tdiff, axis=0) / 2.0 grad = dot(self.inv_cov, diff) result[:,i] = sum(grad * exp(-energy), axis=1)

For some reason I needed to drop the -1 when calculating grad to obtain results coherent with a evaluating the density estimation at p and p+delta in all four directions, which is a sign I might of course be way off here.


  • Jquery UI Sortable, move item automatically
  • cell spacing in div table
  • apply a javascript function to draggable copy
  • Aptana 3 remove bundle (jquery)
  • XSLT foreach repeating nodes to flat
  • How to create a 2D image by rotating 1D vector of numbers around its center element?
  • import scipy.sparse failed
  • netsh acl setting (need alternative method - registry settings?)
  • OSX - always hide certain files
  • Thread 1: EXC_BAD_ACCESS (code =1 address = 0x0)
  • During installation of Django, why do I keep getting ImportError: No module named django?
  • where do I find the xml.dom python package for the python-2.6.0-8.9.28 and I have a suse/x86_64 vers
  • Python pickle not one-to-one: different pickles give same object
  • as3-flash: any way to access all the instances placed in different frames from document class?
  • Tamper-proof configuration files in .NET?
  • Django simple Captcha “No module named fields” error
  • How to use carriage return with multiple line?
  • Combining SpatialPolygonsDataFrame of two neighbour countries
  • Cancel a live stream “fast motion” catch-up in Flash
  • QLineEdit password safety
  • C# - Serializing and deserializing static member
  • Counter field in MS Access, how to generate?
  • Bug in WPF DataGrid
  • Incrementing object id automatically JS constructor (static method and variable)
  • Does CUDA 5 support STL or THRUST inside the device code?
  • WinForms: two way TextBox problem
  • Trying to switch camera back to front but getting exception
  • Javascript + PHP Encryption with pidCrypt
  • Websockets service method fails during R startup
  • Why winpcap requires both .lib and .dll to run?
  • Transpose CSV data with awk (pivot transformation)
  • Free memory of cv::Mat loaded using FileStorage API
  • Why can't I rebase on to an ancestor of source changesets if on a different branch?
  • Angular 2 constructor injection vs direct access
  • how does django model after text[] in postgresql [duplicate]
  • Append folder name and increment by 1 using batch script
  • Programmatically clearing map cache
  • reshape alternating columns in less time and using less memory
  • Android Heatmap on canvas or ImageView
  • Conditional In-Line CSS for IE and Others?