79176

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?

Answer1:

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.

Recommend

  • 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?