Skip to main content
`Ctrl`+`K`
scipy.stats.

# gaussian_kde

class scipy.stats.gaussian_kde(dataset, bw_method=None, weights=None)

  Representation of a kernel-density estimate using Gaussian kernels.

  Kernel density estimation is a way to estimate the probability density function (PDF) of a random variable in a non-parametric way. `gaussian_kde` works for both uni-variate and multi-variate data. It includes automatic bandwidth determination. The estimation works best for a unimodal distribution; bimodal or multi-modal distributions tend to be oversmoothed.

  Parameters:
    dataset array_like

      Datapoints to estimate from. In case of univariate data this is a 1-D array, otherwise a 2-D array with shape (# of dims, # of data).

    bw_method str, scalar or callable, optional

      The method used to calculate the bandwidth factor. This can be ‘scott’, ‘silverman’, a scalar constant or a callable. If a scalar, this will be used directly as factor. If a callable, it should take a `gaussian_kde` instance as only parameter and return a scalar. If None (default), ‘scott’ is used. See Notes for more details.

    weights array_like, optional

      weights of datapoints. This must be the same shape as dataset. If None (default), the samples are assumed to be equally weighted

  Attributes:
    dataset ndarray

      The dataset with which `gaussian_kde` was initialized.

    d int

      Number of dimensions.

    n int

      Number of datapoints.

    neff int

      Effective number of datapoints.

      Added in version 1.2.0.

    factor float

      The bandwidth factor obtained from `covariance_factor`.

    covariance ndarray

      The kernel covariance matrix; this is the data covariance matrix multiplied by the square of the bandwidth factor, e.g. `np.cov(dataset) * factor**2`.

    inv_cov ndarray

      The inverse of covariance.

  Methods

  `evaluate`(points) | Evaluate the estimated pdf on a set of points.
  `__call__`(points) | Evaluate the estimated pdf on a set of points.
  `integrate_gaussian`(mean, cov) | Multiply estimated density by a multivariate Gaussian and integrate over the whole space.
  `integrate_box_1d`(low, high) | Computes the integral of a 1D pdf between two bounds.
  `integrate_box`(low_bounds, high_bounds[, ...]) | Computes the integral of a pdf over a rectangular interval.
  `integrate_kde`(other) | Computes the integral of the product of this kernel density estimate with another.
  `pdf`(x) | Evaluate the estimated pdf on a provided set of points.
  `logpdf`(x) | Evaluate the log of the estimated pdf on a provided set of points.
  `resample`([size, seed]) | Randomly sample a dataset from the estimated pdf.
  `set_bandwidth`([bw_method]) | Compute the bandwidth factor with given method.
  `covariance_factor`() | Computes the bandwidth factor factor.
  `marginal`(dimensions) | Return a marginal KDE distribution

  Notes

  Bandwidth selection strongly influences the estimate obtained from the KDE (much more so than the actual shape of the kernel). Bandwidth selection can be done by a “rule of thumb”, by cross-validation, by “plug-in methods” or by other means; see [3], [4] for reviews. `gaussian_kde` uses a rule of thumb, the default is Scott’s Rule.

  Scott’s Rule [1], implemented as `scotts_factor`, is:

      n**(-1./(d+4)),

  with `n` the number of data points and `d` the number of dimensions. In the case of unequally weighted points, `scotts_factor` becomes:

      neff**(-1./(d+4)),

  with `neff` the effective number of datapoints. Silverman’s suggestion for multivariate data [2], implemented as `silverman_factor`, is:

      (n * (d + 2) / 4.)**(-1. / (d + 4)).

  or in the case of unequally weighted points:

      (neff * (d + 2) / 4.)**(-1. / (d + 4)).

  Note that this is not the same as “Silverman’s rule of thumb” [6], which may be more robust in the univariate case; see documentation of the `set_bandwidth` method for implementing a custom bandwidth rule.

  Good general descriptions of kernel density estimation can be found in [1] and [2], the mathematics for this multi-dimensional implementation can be found in [1].

  With a set of weighted samples, the effective number of datapoints `neff` is defined by:

      neff = sum(weights)^2 / sum(weights^2)

  as detailed in [5].

  `gaussian_kde` does not currently support data that lies in a lower-dimensional subspace of the space in which it is expressed. For such data, consider performing principal component analysis / dimensionality reduction and using `gaussian_kde` with the transformed data.

  References

  [1] (1,2,3)

  D.W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization”, John Wiley & Sons, New York, Chicester, 1992.

  [2] (1,2)

  B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Vol. 26, Monographs on Statistics and Applied Probability, Chapman and Hall, London, 1986.

  [3]

  B.A. Turlach, “Bandwidth Selection in Kernel Density Estimation: A Review”, CORE and Institut de Statistique, Vol. 19, pp. 1-33, 1993.

  [4]

  D.M. Bashtannyk and R.J. Hyndman, “Bandwidth selection for kernel conditional density estimation”, Computational Statistics & Data Analysis, Vol. 36, pp. 279-298, 2001.

  [5]

  Gray P. G., 1969, Journal of the Royal Statistical Society. Series A (General), 132, 272

  [6]

  Kernel density estimation. Wikipedia. https://en.wikipedia.org/wiki/Kernel_density_estimation

  Examples

  Generate some random two-dimensional data:

      >>> import numpy as np
      >>> from scipy import stats
      >>> def measure(n):
      ...     "Measurement model, return two coupled measurements."
      ...     m1 = np.random.normal(size=n)
      ...     m2 = np.random.normal(scale=0.5, size=n)
      ...     return m1+m2, m1-m2

      >>> m1, m2 = measure(2000)
      >>> xmin = m1.min()
      >>> xmax = m1.max()
      >>> ymin = m2.min()
      >>> ymax = m2.max()

  Perform a kernel density estimate on the data:

      >>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
      >>> positions = np.vstack([X.ravel(), Y.ravel()])
      >>> values = np.vstack([m1, m2])
      >>> kernel = stats.gaussian_kde(values)
      >>> Z = np.reshape(kernel(positions).T, X.shape)

  Plot the results:

      >>> import matplotlib.pyplot as plt
      >>> fig, ax = plt.subplots()
      >>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
      ...           extent=[xmin, xmax, ymin, ymax])
      >>> ax.plot(m1, m2, 'k.', markersize=2)
      >>> ax.set_xlim([xmin, xmax])
      >>> ax.set_ylim([ymin, ymax])
      >>> plt.show()

  Compare against manual KDE at a point:

      >>> point = [1, 2]
      >>> mean = values.T
      >>> cov = kernel.factor**2 * np.cov(values)
      >>> X = stats.multivariate_normal(cov=cov)
      >>> res = kernel.pdf(point)
      >>> ref = X.pdf(point - mean).sum() / len(mean)
      >>> np.allclose(res, ref)
      True