Merge pull request #133 from dseuss/feature-custom-meanfn

Feature custom meanfn

Author: cgranade

src/qinfer/clustering.py

@@ -5,7 +5,7 @@
 ##
 # © 2013 Chris Ferrie (csferrie@gmail.com) and
 #        Christopher E. Granade (cgranade@gmail.com)
-#     
+#
 # This file is a part of the Qinfer project.
 # Licensed under the AGPL version 3.
 ##
@@ -46,7 +46,7 @@
 import numpy as np
 import scipy.linalg as la
 
-from qinfer.utils import outer_product, particle_meanfn, particle_covariance_mtx
+from qinfer.utils import outer_product
 from qinfer._exceptions import ResamplerWarning
 import qinfer.metrics as metrics
 
@@ -85,20 +85,20 @@ def particle_clusters(
     of all particles in that cluster. That is, particle ``i`` is in the cluster
     if ``cluster_particles[i] == True``.
     """
-    
-    
+
+
     if weighted == True and particle_weights is None:
         raise ValueError("Weights must be specified for weighted clustering.")
-        
-    # Allocate new arrays to hold the weights and locations.        
+
+    # Allocate new arrays to hold the weights and locations.
     new_weights = np.empty(particle_weights.shape)
     new_locs    = np.empty(particle_locations.shape)
-    
+
     # Calculate and possibly reweight the metric.
     if weighted:
         M = sklearn.metrics.pairwise.pairwise_distances(particle_locations, metric=metric)
         M = metrics.weighted_pairwise_distances(M, particle_weights, w_pow=w_pow)
-    
+
         # Create and run a SciKit-Learn DBSCAN clusterer.
         clusterer = sklearn.cluster.DBSCAN(
             min_samples=min_particles,
@@ -113,33 +113,33 @@ def particle_clusters(
             metric=metric
         )
         cluster_labels = clusterer.fit_predict(particle_locations)
-    
+
     # Find out how many clusters were identified.
     # Cluster counting logic from:
     # [http://scikit-learn.org/stable/auto_examples/cluster/plot_dbscan.html].
     is_noise = -1 in cluster_labels
     n_clusters = len(set(cluster_labels)) - (1 if is_noise else 0)
-    
+
     # If more than 10% of the particles were labeled as NOISE,
     # warn.
     n_noise = np.sum(cluster_labels == -1)
     if n_noise / particle_weights.shape[0] >= 0.1:
         warnings.warn("More than 10% of the particles were classified as NOISE. Consider increasing the neighborhood size ``eps``.", ResamplerWarning)
-    
+
     # Print debugging info.
     if not quiet:
         print("[Clustering] DBSCAN identified {} cluster{}. "\
               "{} particles identified as NOISE.".format(
                   n_clusters, "s" if n_clusters > 1 else "", n_noise
               ))
-        
+
     # Loop over clusters, calling the secondary resampler for each.
     # The loop should include -1 if noise was found.
     for idx_cluster in range(-1 if is_noise else 0, n_clusters):
         # Grab a boolean array identifying the particles in a  particular
         # cluster.
         this_cluster = cluster_labels == idx_cluster
-        
+
         yield idx_cluster, this_cluster
-    
+
 

src/qinfer/distributions.py

@@ -45,6 +45,7 @@
 from qinfer import utils as u
 from qinfer.metrics import rescaled_distance_mtx
 from qinfer.clustering import particle_clusters
+from qinfer._exceptions import ApproximationWarning
 
 import warnings
 
@@ -321,6 +322,70 @@ def sample(self, n=1):
         ), len(cumsum_weights) - 1)]
 
     ## MOMENT FUNCTIONS ##
+    @staticmethod
+    def particle_mean(weights, locations):
+        r"""
+        Returns the arithmetic mean of the `locations` weighted by `weights`
+
+        :param numpy.ndarray weights: Weights of each particle in array of
+            shape ``(n_particles,)``.
+        :param numpy.ndarray locations: Locations of each particle in array
+            of shape ``(n_particles, n_modelparams)``
+        :rtype: :class:`numpy.ndarray`, shape ``(n_modelparams,)``.
+        :returns: An array containing the mean
+        """
+        return np.dot(weights, locations)
+
+    @classmethod
+    def particle_covariance_mtx(cls, weights, locations):
+        """
+        Returns an estimate of the covariance of a distribution
+        represented by a given set of SMC particle.
+
+        :param weights: An array of shape ``(n_particles,)`` containing
+            the weights of each particle.
+        :param location: An array of shape ``(n_particles, n_modelparams)``
+            containing the locations of each particle.
+        :rtype: :class:`numpy.ndarray`, shape
+            ``(n_modelparams, n_modelparams)``.
+        :returns: An array containing the estimated covariance matrix.
+        """
+        # Find the mean model vector, shape (n_modelparams, ).
+        mu = cls.particle_mean(weights, locations)
+
+        # Transpose the particle locations to have shape
+        # (n_modelparams, n_particles).
+        xs = locations.transpose([1, 0])
+        # Give a shorter name to the particle weights, shape (n_particles, ).
+        ws = weights
+
+        cov = (
+            # This sum is a reduction over the particle index, chosen to be
+            # axis=2. Thus, the sum represents an expectation value over the
+            # outer product $x . x^T$.
+            #
+            # All three factors have the particle index as the rightmost
+            # index, axis=2. Using the Einstein summation convention (ESC),
+            # we can reduce over the particle index easily while leaving
+            # the model parameter index to vary between the two factors
+            # of xs.
+            #
+            # This corresponds to evaluating A_{m,n} = w_{i} x_{m,i} x_{n,i}
+            # using the ESC, where A_{m,n} is the temporary array created.
+            np.einsum('i,mi,ni', ws, xs, xs)
+            # We finish by subracting from the above expectation value
+            # the outer product $mu . mu^T$.
+            - np.dot(mu[..., np.newaxis], mu[np.newaxis, ...])
+        )
+
+        # The SMC approximation is not guaranteed to produce a
+        # positive-semidefinite covariance matrix. If a negative eigenvalue
+        # is produced, we should warn the caller of this.
+        assert np.all(np.isfinite(cov))
+        if not np.all(la.eig(cov)[0] >= 0):
+            warnings.warn('Numerical error in covariance estimation causing positive semidefinite violation.', ApproximationWarning)
+
+        return cov
 
     def est_mean(self):
         """
@@ -329,13 +394,8 @@ def est_mean(self):
         :rtype: :class:`numpy.ndarray`, shape ``(n_mps,)``.
         :returns: An array containing the an estimate of the mean model vector.
         """
-        return np.sum(
-            # We need the particle index to be the rightmost index, so that
-            # the two arrays align on the particle index as opposed to the
-            # modelparam index.
-            self.particle_weights * self.particle_locations.transpose([1, 0]),
-            axis=1
-        )
+        return self.particle_mean(self.particle_weights,
+                                  self.particle_locations)
 
     def est_meanfn(self, fn):
         """
@@ -372,9 +432,8 @@ def est_covariance_mtx(self, corr=False):
         :returns: An array containing the estimated covariance matrix.
         """
 
-        cov = u.particle_covariance_mtx(
-            self.particle_weights,
-            self.particle_locations)
+        cov = self.particle_covariance_mtx(self.particle_weights,
+                                           self.particle_locations)
 
         if corr:
             dstd = np.sqrt(np.diag(cov))
@@ -448,8 +507,8 @@ def est_cluster_moments(self, cluster_opts=None):
             yield (
                 cluster_label,
                 sum(w), # The zeroth moment is very useful here!
-                u.particle_meanfn(w, l, lambda x: x),
-                u.particle_covariance_mtx(w, l)
+                self.particle_mean(w, l),
+                self.particle_covariance_mtx(w, l)
             )
 
     def est_cluster_covs(self, cluster_opts=None):
@@ -467,7 +526,7 @@ def est_cluster_covs(self, cluster_opts=None):
         ws = cluster_moments['weight'][:, np.newaxis, np.newaxis]
 
         within_cluster_var = np.sum(ws * cluster_moments['cov'], axis=0)
-        between_cluster_var = u.particle_covariance_mtx(
+        between_cluster_var = self.particle_covariance_mtx(
             # Treat the cluster means as a new very small particle cloud.
             cluster_moments['weight'], cluster_moments['mean']
         )