Bandwidth parameters in multivariate KDE using scipy.stats.gaussian_kde

Question

I am working on a project which involves implementing in Python two different density estimation functions over multivariate data; one using N-d histograms and the other using kernel density estimation (KDE).

I used scipy.stats.gaussian_kde, and realised that as the dimensions differed in variance in my underlying data, the KDE function was less able to accurately estimate the underlying probability density.

It seems that we need to have a different bandwidth parameter for each dimension (manifest, I gather, as a bandwidth matrix). scipy.stats.gaussian_kde does not appear to support bandwidth matrix arguments, so my questions are:

1. Am I right in thinking that this function is flawed/useless without performing some normalisation-type operation on the data to ensure SD is somehow the same for each dimension? If so, it seems like a glaring problem in the implementation.

2. What is the simplest Python alternative to plug into my code, to calculate KDE given a bandwidth matrix? (Ideally with a similar interface to avoid much refactoring)

Answer 1

Multivariate kernel density estimation can be defined in terms of

1. product of univariate kernels $K^P(\mathrm{x}) = \prod_{i=1}^{d} \kappa(x_i)$,
2. as a symmetric kernel $K^S(\mathrm{x}) \propto\kappa\{(\mathrm{x}'\mathrm{x})^{1/2}\}$,
3. or in terms of standalone multivariate kernel, e.g. multivariate Gaussian distribution.

There are also different possible choices of bandwidth matrix,

1. it can have equal bandwidth for each of the variables $\mathrm{H} = h^2\mathrm{I}_d$,
2. different for different variables $\mathrm{H} = \mathrm{diag}(h_1^2, h_2^2, \dots, h_d^2)$,
3. or it could be a covariance matrix.

The three choices are illustrated by Wand and Jones in their Kernel Smoothing book using the following figure of two-dimensional case.

The first choice is "symmetric", it assumes no correlation and equal variances. Second allows for unequal variances. Third allows additionally for correlation between variables.

The scipy documentation does not tell us much about the kind of multivariate kernel that they are using. It only tells us that it uses Scotts or Silvermans rules of thumb for selecting the bandwidth, so it estimates some constant $h^2$ and the covariance matrix is either same for all variables or is a scaling factor for covariance matrix (more likely, but you'd need to check the source code). Nonetheless, scipy is using rule of a thumb for choosing the bandwidth, so this does not have to be optimal choice and I'd encourage you to look for packages that implement more sophisticted approaches (R has several, I can't tell for python).

Answer 2

You can provide your own covariance matrix by inheriting from the gaussian_kde class as follows:

from scipy.stats._kde import gaussian_kde
import numpy as np

class custom_kde(gaussian_kde):
    def __init__(self, dataset, bw_method=None, weights=None, covariance=None):
        self.covariance = covariance
        super().__init__(dataset, bw_method=bw_method, weights=weights)        
    
    def _compute_covariance(self):    
        # Copied from the base gaussian_kde class, except for the covariance part    
        self.factor = self.covariance_factor()        
        if not hasattr(self, '_data_inv_cov'):
            self._data_covariance = np.atleast_2d(np.cov(self.dataset, rowvar=1, bias=False, aweights=self.weights))
            self._data_inv_cov = np.linalg.inv(self._data_covariance)
    
        if self.covariance is None:
            self.covariance = self._data_covariance * self.factor**2
        self.inv_cov = np.linalg.inv(self.covariance)
        L = np.linalg.cholesky(self.covariance*2*np.pi)
        self.log_det = 2*np.log(np.diag(L)).sum()

Example usage:

kernel = custom_kde(X_arr, Y_arr, covariance=np.array([[1.0, 0.5], [0.5, 1.0]]))

Answer 3

I was looking for an answer to this bandwidth matrix optimisation problem and I found this excellent other thread, so I thought I'd drop it here:

Difference in bandwidth for scikit-learn KDE and multivariate KDE of statsmodels.

In short, it says that the sklearn KernelDensity() implementation uses bandwidth as a multiplier of the diagonal matrix (so second case of Tim's answer), while statsmodel's KDEMultivariate() estimates different multipliers (so third picture, I believe). I am not sure how this compares to scipy which multiplies the covariance matrix by the single scalar. It looks to fall in the same case as KDEMultivariate(), but with a little less control over the dimension-specific toggling. From what I understand (again from that other stackoverflow answer), they both use rule of thumb for coming up with the covariance matrix.

Answer 4

This will allow you to manually specify a covariance matrix for your kernel, made by stripping Joe's answer down to the bare bones.

from scipy.stats import gaussian_kde as kde

class custom_kde(kde):
    def __init__(self, dataset, covariance):
        self.covariance = covariance
        super().__init__(dataset, bw_method=1.0, weights=None)

    def _compute_covariance(self):
        self.inv_cov = np.linalg.inv(self.covariance)

Example use

from scipy.stats import multivariate_normal as mvn
from matplotlib import pyplot as plt
plt.ion()

# covariance matrix with std devs of 10 and 1 in the x and y
# and a correlation coefficient of 5/(10*1) = 0.5
covariance = np.array([[100, 5], [5, 1]])

# compute KDE on a grid based on a single data point 
# with the manually set covariance
data = [[0], [0]]
kernel = custom_kde(data, covariance)
x, y = np.mgrid[-30:30:100j, -3:3:100j]
xy = np.vstack([x.ravel(), y.ravel()])
density = kernel(xy)
density = density.reshape(x.shape)

# plot it up
plt.imshow(np.rot90(density), cmap=plt.get_cmap('YlOrRd'),
              extent=(-30, 30, -3, 3))
plt.gca().set_aspect('auto')

A multivariate Gaussian with the covariances we specified is what we get.