# Bandwidth parameters in multivariate KDE using scipy.stats.gaussian_kde

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)

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).

• Thanks - I have been passing in a scalar bandwidth parameter to scipy's gaussian_kde. Your answer helped me spot that this is in fact applied as an element-wise multiplier to a covariance bandwidth matrix - corresponding to your third option. (Docs for scipy.stats.gaussian_kde.covariance_factor). However - I was hoping to find a way to define this matrix in order to deal 'optimally' with multidimensional data with different variances in each dimension. My research suggests only statsmodels.nonparametric.kernel_density.KDEMultivariate provides this, but I haven't worked out yet how to use it!
– Zac
Jul 20, 2018 at 22:25

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]]))