U.S. flag

An official website of the United States government
bit Xplor-NIH VMD-XPLOR
Pydoc: module densityEstimation
 
densityEstimation		 index


 
Module with utilities for density estimation.
 
When using the kernel density estimation capabilities of this module, please
cite:
 
Bermejo, G.A., Clore, G.M., and Schwieters, C.D. (2012). Smooth statistical
torsion angle potential derived from a large conformational database via
adaptive kernel density estimation improves the quality of NMR protein
structures. Protein Sci. 21, 1824-1836.

 
Classes
       
builtins.object
KDE
NormalKDE

 
class KDE(builtins.object)
    KDE(dataset, width=None, localwidths=None, symmetric=True, cov=None, wrapped=False)
 
General kernel density estimate with a fixed- or variable-width kernel.
 
Subclass and implement the method evaluate according to a specific kernel
function (e.g., class NormalKDE(KDE) evaluates the estimates using a
multivariate normal kernel).  Optionally, the kernel-specific method
optimal_width, that calculates an optimal width assuming an underlying
probability density function (pdf), can also be implemented.
 
The general formula used for density estimation is that given by (Wand and
Jones, 1993), modified to include local widths as in (Silverman, 1986):
 
f(x) = n^-1 SUM_i{h^-d |C|^-0.5 lambdai^-d K[(h lambdai)^-1 C^-0.5 (x-Xi)]}
 
where x and Xi are d-dimensional vectors (Xi represents sample data point
i), n is the total number of points, C is a symmetric positive definite
d-by-d matrix with determinant |C|, and h and lambdai are, respectively, the
global and local widths of kernel K at Xi.  The following discussion will
assume kernel K is a symmetric pdf (e.g., multivariate standard normal), and
refer to the expression within the sum as the "effective kernel" (which
might not be symmetric).
 
The constructor requires a dataset, a cdsMatrix.CDSMatrix_type(d, n)
(type=int, double) that contains the sample points from which the density is
estimated.  For example, for two 3D points with coordinates (1.0, 2.0, 3.0)
and (4.0, 5.0, 6.0):
 
    dataset = cdsMatrix.CDSMatrix_double(3, 2).fromList([[1.0, 4.0],
                                                         [2.0, 5.0],
                                                         [3.0, 6.0]])
 
The constructor call is something like:
 
    kde = NormalKDE(dataset)
    
Optionally, the arguments width, localwidths, symmetric, cov, and wrapped
can also be specified.  width is the global kernel width h; if the method
optimal_width is implemented, an optimal width is automatically calculated
when width is not supplied.  localwidths is a sequence of length n that
specifies a local width for each kernel used in the density estimation:
localwidths[i]=lambdai; if not supplied lambdai=1 for all i (i.e., a fixed-
width effective kernel is used).  Different combinations of values for
arguments symmetric and cov determine matrix C:
 
symmetric=True => C is the identity matrix;
    
symmetric=False and cov=None => C is the covariance matrix of the data
(the Fukunaga method; see Silverman 1986, Sec. 4.2.1, p. 77);
    
symmetric=False and cov=A => C is the arbitrary symmetric positive definite
matrix A, input as either a cdsMatrix.CDSMatrix_type(d, d) (type=int,
double) or a sequence of diagonal elements for a diagonal A.
 
Thus, symmetric=True yields a radially symmetric effective kernel, while
options with symmetric=False yield an asymmetric one (except for trivial
cases such as setting cov=A=constant*identity matrix). 
The remaining argument, wrapped, specifies whether the input dataset has
been augmented by surrounding it with one copy of it (e.g., by using
function wraparound in this module).  
 
The evaluate method in a KDE subclass should be coded so that density
estimates at arbitrary points are obtained via:
 
    vals = kde.evaluate(pts)
    
pts is a cdsMatrix.CDSMatrix_type(d, m) (type=int, double), where m is
the total number of points.  For example, if kde is to be evaluated at
(1.5, 2.5, 3.5), (4.5, 5.5, 6.5), and (7.5, 8.5, 9.5):
 
    pts = cdsMatrix.CDSMatrix_double(3, 3).fromList([[1.5, 4.5, 7.5],
                                                     [2,5, 5.5, 8.5],
                                                     [3.5, 6.5, 9.5]])
                                                     
The returned values, vals, is a cdsVector.CDSVector_double of size m,
where vals[i] is the density at the point in column i of pts (e.g., vals[1]
is the density at (4.5, 5.5, 6.5) in the above example).
 
References:
Silverman, B.W. (1986).  Density Estimation for Statistics and Data
Analysis.  Chapman & Hall.
Wand, M.P. and Jones, M.C. (1993).  Comparison of Smoothing Parametrizations
in Bivariate Kernel Density-Estimation. J. Amer. Statistical Assoc. 88(422):
520-528.
 
  Methods defined here:
__init__(self, dataset, width=None, localwidths=None, symmetric=True, cov=None, wrapped=False)
Initialize self.  See help(type(self)) for accurate signature.
evaluate(self, points)
Return kernel-specific density estimates at arbitrary points.
get_samplecov(self)
Return a cdsMatrix.CDSMatrix_double with the sample covariance matrix.
get_widths(self)
Return a cdsVector.CDSVector_double of localwidth*width at each point.
ndim(self)
Return the dimensionality of the estimated density function.
npts(self)
Return the number of sample points used for density estimation.
optimal_width(self)
Return kernel-specific optimal width assumming an underlying pdf.
 
The implementaton of this method in a subclass is not needed if width is
to be set only subjectively (i.e., specified upon intance creation).
Raise ValueError if this method is called while not set in the subclass.
Data descriptors defined here:
__dict__

 
dictionary for instance variables (if defined)
__weakref__

 
list of weak references to the object (if defined)

 
class NormalKDE(KDE)
    NormalKDE(dataset, width=None, localwidths=None, symmetric=True, cov=None, wrapped=False)
 
Kernel density estimate with a multivariate normal kernel.
 
It uses kernel
 
    K(x) = (2 pi)^-(d/2) exp(-0.5 x^T x)
 
where x is a d-dimensional vector with transpose x^T.  It estimates the
density using the general formula given in the KDE class.  See KDE class for
a detailed description of the interface.
 
 
Method resolution order:
NormalKDE
KDE
builtins.object
Methods defined here:
evaluate(self, points)
Return a cdsVector.CDSVector_double of density estimates at points.
 
points is a cdsMatrix.CDSMatrix_type(d, m) (type=int, double), where
m is the total number of points.  For example, a NormalKDE instance is
evaluated at (1.5, 2.5, 3.5), (4.5, 5.5, 6.5), and (7.5, 8.5, 9.5) as
follows:
 
points = cdsMatrix.CDSMatrix_double(3, 3).fromList([[1.5, 4.5, 7.5],
                                                    [2.5, 5.5, 8.5],
                                                    [3.5, 6.5, 9.5]])
                                                     
The returned values vector, vals, has size m, and vals[i] is the density
estimate at the point in column i of the input points matrix (e.g.,
vals[1] is the density at (4.5, 5.5, 6.5) in the above example).
marginalize(self, dims)
Return a NormalKDE instance with the marginal density over dims.
 
dims argument is either an integer or a sequence of intergers specifying
the dimension(s) over which to marginalize.
optimal_width(self)
Return optimal width assuming the underlying pdf is normal.
 
If symmetric=False it returns hopt = [4/n(2d+1)]^[1/(d+4)].
If symmetric=True it returns sqrt(Tr{S}/d)*hopt, where Tr{S} is the
trace of the covariance matrix of the data (Silverman, 1986, Sec. 4.3.2,
pp.86, 87).
 
Reference: Silverman, B.W. (1986). Density Estimation for Statistics
and Data Analysis. Chapman & Hall.
Methods inherited from KDE:
__init__(self, dataset, width=None, localwidths=None, symmetric=True, cov=None, wrapped=False)
Initialize self.  See help(type(self)) for accurate signature.
get_samplecov(self)
Return a cdsMatrix.CDSMatrix_double with the sample covariance matrix.
get_widths(self)
Return a cdsVector.CDSVector_double of localwidth*width at each point.
ndim(self)
Return the dimensionality of the estimated density function.
npts(self)
Return the number of sample points used for density estimation.
Data descriptors inherited from KDE:
__dict__

 
dictionary for instance variables (if defined)
__weakref__

 
list of weak references to the object (if defined)

 
Functions
       
adaptive_kde(dataset, width=None, symmetric=True, cov=None, wrapped=False, alpha=0.5, pilot=<class 'densityEstimation.NormalKDE'>, adaptive=<class 'densityEstimation.NormalKDE'>)
Return a class 'adaptive' instance: the adaptive kernel density estimate.
 
Arguments dataset, width, symmetric, cov, and wrapped are those required for
KDE instantiation; alpha is the sensitivity parameter used in the local
width calculation from the pilot density estimate:
 
    lambdai = [pilot_f(Xi)/g]^-alpha
 
with lambdai the local width of the kernel centered at vector Xi, where the
pilot desity is pilot_f(Xi), and log[g] = n^-1 SUM_i{log[pilot_f(Xi)]} (eq.
5.7 in Silverman 1986, Sec. 5.3, p.101).  The kernels used in the pilot and
adaptive stages are respectively specified with arguments pilot and
adaptive, which expect the names of the corresponding classes; the global
width of both kernels is the same, specified with the width argument.  The
overall procedure is described by (Silverman, 1986, Sec. 5.3); for more
details see KDE docstring.
 
Reference: Silverman, B.W. (1986).  Density Estimation for Statistics and
Data Analysis.  Chapman & Hall.
gengrid(npoints, ranges, periodic=False)
Return a Python itertools.product object representing a d-dimesional grid.
 
Note: This is a Python generator version of the grid function; another
difference with grid is that it lacks the functionality associated with the
axes argument (there is no such argument here).
 
npoints is a sequence with the number of points in each dimension.  ranges
is a sequence of (min, max) pairs, with the min and max values in each
dimension (the order of values within a pair is irrelevant).  The number of
dimensions, d, is len(npoints) (or len(ranges); ValueError is raised if
len(npoints)!=len(ranges)).  For example:
 
    npoints = (5, 3)
    ranges = ([-10, 10], [-5, 5])
 
represents a system with d=2, where the first dimension or axis has values
in the [-10, 10] range, and the second axis in the [-5, 5].  5 points are
requested from the first axis and 3 from the second; points are equally
spaced to make the grid, i.e.
 
    Axis 1: [-10.0, -5.0, 0.0, 5.0, 10.0]
    Axis 2: [-5.0, 0.0, 5.0]
 
The grid associated with the above axes is the matrix
 
      -10, -10, -10, -5, -5, ..., 10      # from axis 1
       -5,   0,   5, -5,  0, ...,  5      # from axis 2
 
The matrix columns are the multidimensional (in this example two-
dimensional) points in the grid.  They are returned via a Python
itertools.product object, a generator that yields each point packed into a
tuple.
                                                    
The optional argument periodic is a sequence with True/False for each axis;
if True, the last point in the corresponding axis is omitted.  For example,
gengrid(npoints, ranges, (True, False)) is associated with a grid generated
from axes
 
    Axis 1: [-10.0, -5.0, 0.0, 5.0]
    Axis 2: [-5.0, 0.0, 5.0]
 
Alternatively, all axes can be given the same True or False periodic value
by setting periodic=True or periodic=False, respectively, i.e.,
 
gengrid(npoints, ranges, (True, True)) <-> gengrid(npoints, ranges, True) 
gengrid(npoints, ranges, (False, False)) <-> gengrid(npoints, ranges, False)
grid(npoints, ranges, periodic=False, axes=False)
Return a cdsMatrix.CDSMatrix_double representing a d-dimesional grid.
 
npoints is a sequence with the number of points in each dimension.  ranges
is a sequence of (min, max) pairs, with the min and max values in each
dimension (the order of values within a pair is irrelevant).  The number of
dimensions, d, is len(npoints) (or len(ranges); ValueError is raised if
len(npoints)!=len(ranges)).  For example:
 
    npoints = (5, 3)
    ranges = ([-10, 10], [-5, 5])
 
represents a system with d=2, where the first dimension or axis has values
in the [-10, 10] range, and the second axis in the [-5, 5].  5 points are
requested from the first axis and 3 from the second; points are equally
spaced to make the grid, i.e.
 
    Axis 1: [-10.0, -5.0, 0.0, 5.0, 10.0]
    Axis 2: [-5.0, 0.0, 5.0]
 
The grid can be generated via grid(npoints, ranges), which, following the
above example, returns a cdsMatrix.CDSMatrix_double(2, 5*3) representing
 
    [[-10, -10, -10, -5, -5, ..., 10]   # axis 1
     [ -5,   0,   5, -5,  0, ...,  5]]  # axis 2
                                                    
The optional argument periodic is a sequence with True/False for each axis;
if True, the last point in the corresponding axis is omitted.  For example,
grid(npoints, ranges, (True, False)) returns a grid with 
 
    Axis 1: [-10.0, -5.0, 0.0, 5.0]
    Axis 2: [-5.0, 0.0, 5.0]
 
Alternatively, all axes can be given the same True or False periodic value
by setting periodic=True or periodic=False, respectively, i.e.,
 
    grid(npoints, ranges, (True, True)) <-> grid(npoints, ranges, True) 
    grid(npoints, ranges, (False, False)) <-> grid(npoints, ranges, False)
 
The axes may be obtained by setting argument axes=True, in which case the
function returns a tuple with the grid (the cdsMatrix.CDSMatrix_double)
as first element and the axes, packed into a list, as second element.  For
example, grid(npoints, ranges, True, True) returns
 
    (grid(npoints, ranges, True), [[-10.0, -5.0, 0.0, 5.0], [-5.0, 0.0]])
 
i.e., (the grid, [axis 1, axis 2]).
wraparound(dataset, spectral_widths)
Return an augmented dataset by adding shifted copies of it.
 
Impose a periodic or "wrap around" boundary condition on the dataset by
augmenting it via the addition of shifted copies of it (Silverman, 1986,
Sec. 2.10, p. 31, 32).
 
dataset is cdsMatrix.CDSMatrix_type(d, n) (type=int, double), where d is
the dimensionality and n the number of points (like in KDE instantiation).
For example, two 3D points with coordinates (1.0, 2.0, 3.0) and
(4.0, 5.0, 6.0) yield:
 
    dataset = cdsMatrix.CDSMatrix_double(3, 2).fromList([[1.0, 4.0],
                                                         [2.0, 5.0],
                                                         [3.0, 6.0]])
                                                         
After wrapping around the number of points is (3^d)n.
spectral_widths is a sequence with the size of the domain in each dimension.
For example, if dataset contains directions or angles defined within the
range [-pi, pi], the associated size in spectral_widths is 2pi.  Sizes in
spectral_widths and dimensions in dataset are matched by offset: row i in
dataset corresponds to spectral_widths[i]; ValueError is raised if the
dimensionality of dataset and spectral_widths do not match.
 
Reference: Silverman, B.W. (1986).  Density Estimation for Statistics and
Data Analysis.  Chapman & Hall.