Gaussian Mixture Models (#285)

Author: jonathf

CHANGELOG.rst

@@ -1,6 +1,10 @@
 Master Branch
 =============
 
+ADDED:
+  * Gaussian Mixture Model: `GaussianMixture`.
+  * Tutorial for how to use `scikit-learn` mixture models to fit a model, and
+    `chaospy` to generate quasi-random samples and orthogonal polynomials.
 CHANGED:
   * `chaospy.Trunc` updated to take both `lower` and `upper` at the same time.
 

chaospy/distributions/__init__.py

@@ -73,7 +73,7 @@
 from .operators import *
 from .constructor import construct
 from .approximation import *
-from .kernel import GaussianKDE
+from .kernel import GaussianKDE, GaussianMixture
 
 from . import (
     baseclass, sampler, approximation,

chaospy/distributions/kernel/__init__.py

@@ -1,2 +1,3 @@
 """Kernel density estimation."""
 from .gaussian import GaussianKDE
+from .mixture import GaussianMixture

chaospy/distributions/kernel/baseclass.py

@@ -78,8 +78,8 @@ def __init__(
             dependency_type="accumulate", length=len(samples),
         )
 
-        self.samples = samples
-        self.h_mat = covariance
+        self._samples = samples
+        self._covariance = covariance
         self._permute = numpy.eye(len(rotation), dtype=int)[rotation]
         self._pcovariance = numpy.matmul(numpy.matmul(
             self._permute, covariance), self._permute.T)
@@ -107,9 +107,9 @@ def get_parameters(self, idx, cache, assert_numerical=True):
 
     def _pdf(self, x_loc, idx, dim, cache):
         """Kernel density function."""
-        s, t = numpy.mgrid[:x_loc.shape[-1], :self.samples.shape[-1]]
+        s, t = numpy.mgrid[:x_loc.shape[-1], :self._samples.shape[-1]]
         if not dim:
-            samples = self.samples[idx, t]
+            samples = self._samples[idx, t]
             z_loc = ((x_loc[s]-samples)*self._fwd_transform[:, 0, 0])
             self._zloc = z_loc[:, :, numpy.newaxis]
             kernel = self._kernel(self._zloc)/self._inv_transform[:, 0, 0]
@@ -122,7 +122,7 @@ def _pdf(self, x_loc, idx, dim, cache):
             x_loc = [self._get_cache(dim_, cache, get=0)
                      for dim_ in self._rotation[:dim]] + [x_loc]
             x_loc = numpy.dstack([x[s] for x in x_loc])
-            samples = numpy.dstack([self.samples[dim_, t]
+            samples = numpy.dstack([self._samples[dim_, t]
                                     for dim_ in self._rotation[:dim+1]])
             zloc = numpy.sum((x_loc-samples)*self._fwd_transform[:, dim, :dim+1], -1)
             self._zloc = numpy.dstack([self._zloc[:, :, :dim], zloc])
@@ -138,9 +138,9 @@ def _pdf(self, x_loc, idx, dim, cache):
 
     def _cdf(self, x_loc, idx, dim, cache):
         """Forward mapping."""
-        s, t = numpy.mgrid[:x_loc.shape[-1], :self.samples.shape[-1]]
+        s, t = numpy.mgrid[:x_loc.shape[-1], :self._samples.shape[-1]]
         if not dim:
-            z_loc = (x_loc[s]-self.samples[idx, t])*self._fwd_transform[:, 0, 0]
+            z_loc = (x_loc[s]-self._samples[idx, t])*self._fwd_transform[:, 0, 0]
             self._zloc = z_loc[:, :, numpy.newaxis]
             out = numpy.sum(self._ikernel(z_loc)*self.weights, axis=-1)
             assert out.shape == x_loc.shape, (out.shape, x_loc.shape)
@@ -150,7 +150,7 @@ def _cdf(self, x_loc, idx, dim, cache):
                      for dim_ in self._rotation[:dim]] + [x_loc]
             x_loc = numpy.dstack([x[s] for x in x_loc])
 
-            samples = numpy.dstack([self.samples[dim_, t]
+            samples = numpy.dstack([self._samples[dim_, t]
                                     for dim_ in self._rotation[:dim+1]])
             zloc = numpy.sum((x_loc-samples)*self._fwd_transform[:, dim, :dim+1], -1)
             self._zloc = numpy.dstack([self._zloc[:, :, :dim], zloc])
@@ -161,7 +161,7 @@ def _cdf(self, x_loc, idx, dim, cache):
 
     def _ppf(self, u_loc, idx, dim, cache):
         """Inverse mapping."""
-        xloc0 = numpy.quantile(self.samples[idx], u_loc)
+        xloc0 = numpy.quantile(self._samples[idx], u_loc)
         out = chaospy.approximate_inverse(
             self, idx, u_loc, xloc0=xloc0, cache=cache, iterations=1000)
         return out

chaospy/distributions/kernel/gaussian.py

@@ -46,6 +46,14 @@ class GaussianKDE(KernelDensityBaseclass):
 
     """
 
+    @property
+    def samples(self):
+        return self._samples
+
+    @property
+    def h_mat(self):
+        return self._covariance
+
     @staticmethod
     def _kernel(z_loc):
         """The underlying density kernel."""

chaospy/distributions/kernel/mixture.py

@@ -0,0 +1,78 @@
+"""Gaussian Mixture Model."""
+import numpy
+
+from .gaussian import GaussianKDE
+
+
+class GaussianMixture(GaussianKDE):
+    """
+    Gaussian Mixture Model.
+
+    A Gaussian mixture model is a probabilistic model that assumes all the data
+    points are generated from a mixture of a finite number of Gaussian
+    distributions with unknown parameters. One can think of mixture models as
+    generalizing K-means clustering to incorporate information about the
+    covariance structure of the data as well as the centers of the latent
+    Gaussians.
+
+    Attributes:
+        means:
+            Sequence of means.
+        covariances:
+            Sequence of covariance matrices.
+        weights:
+            How much each sample is weighted. Either a scalar when the samples
+            are equally weighted, or a vector with the same length as the
+            number of mixed models.
+
+    Examples:
+        >>> means = ([0, 1], [1, 0])
+        >>> covariances = ([[1, 0], [0, 1]], [[1, 0.5], [0.5, 1]])
+        >>> distribution = GaussianMixture(means, covariances)
+        >>> uloc = [[0, 0, 1, 1], [0, 1, 0, 1]]
+        >>> distribution.pdf(uloc).round(4)
+        array([0.0954, 0.092 , 0.1212, 0.0954])
+        >>> distribution.fwd(uloc).round(4)
+        array([[0.3293, 0.3293, 0.6707, 0.6707],
+               [0.3699, 0.6731, 0.3711, 0.734 ]])
+        >>> distribution.inv(uloc).round(4)
+        array([[-8.9681, -8.9681,  8.0521,  8.0521],
+               [-9.862 , 10.1977, -9.5929, 10.2982]])
+        >>> distribution.mom([(0, 1, 1), (1, 0, 1)]).round(4)
+        array([0.5 , 0.5 , 0.25])
+
+    """
+
+    @property
+    def means(self):
+        return self._samples.T
+
+    @property
+    def covariances(self):
+        return numpy.swapaxes(self._covariance, 0, 2)
+
+    def __init__(self, means, covariances, weights=None, rotation=None):
+        """
+        Args:
+            means (numpy.ndarray):
+                Sequence of mean values. With shape `(n_components, n_dim)`.
+            covariances (numpy.ndarray):
+                Sequence of covariance matrices.
+                With shape `(n_components, n_dim, n_dim)`.
+            weights (Optional[numpy.ndarray]):
+                Weights of the samples. This must have the shape
+                `(n_components,)`. If omitted, each sample is assumed to be
+                equally weighted.
+
+        """
+        means = numpy.atleast_2d(numpy.transpose(means))
+        n, m = means.shape
+
+        covariances = numpy.asfarray(covariances)
+        if covariances.ndim in (1, 2):
+            covariances = numpy.broadcast_to(covariances.T, (n, n, m))
+        else:
+            covariances = numpy.swapaxes(covariances, 0, 2)
+        assert covariances.shape == (n, n, m)
+        super(GaussianMixture, self).__init__(
+            samples=means, h_mat=covariances, weights=weights, rotation=rotation)

docs/tutorials/README.rst

@@ -13,6 +13,7 @@ Tutorials
    advanced/scikitlearn_regression
    advanced/intrusive_galerkin
    advanced/kernel_density_estimation
+   advanced/gaussian_mixture_model
    polynomial/expansion_construction
    polynomial/polynomial_evaluation
    polynomial/wiener_askey_scheme
@@ -76,10 +77,14 @@ Advanced Topics
    function of a random variable in a non-parametric way. It works for both
    uni-variate and multi-variate data. It includes automatic bandwidth
    determination.
+`Gaussian Mixture Model`_
+   A Gaussian mixture model is a probabilistic model constructed from a mixture
+   of a finite number of Gaussian distributions.
 
 .. _Scikit-Learn Regression: ./advanced/scikitlearn_regression.ipynb
 .. _Intrusive Galerkin: ./advanced/intrusive_galerkin.ipynb
 .. _Kernel Density Estimation: ./advanced/kernel_density_estimation.ipynb
+.. _Gaussian Mixture Model: ./advanced/gaussian_mixture_model.ipynb
 
 Polynomial Behavior
 -------------------