def est_cluster_moments(self, cluster_opts=None):
        # TODO: document

        if cluster_opts is None:
            cluster_opts = {}

        for cluster_label, cluster_particles in particle_clusters(
                self.particle_locations, self.particle_weights,
                **cluster_opts
            ):

            w = self.particle_weights[cluster_particles]
            l = self.particle_locations[cluster_particles]
            yield (
                cluster_label,
                sum(w), # The zeroth moment is very useful here!
                self.particle_mean(w, l),
                self.particle_covariance_mtx(w, l)
            )