Modified to track sample variance for proper calculation of explained_variance_ratio

kastnerkyle · kastnerkyle · commit e0cb293ef936 · 2014-07-01T12:47:10.000+02:00
diff --git a/sklearn/decomposition/incremental_pca.py b/sklearn/decomposition/incremental_pca.py
@@ -11,12 +11,33 @@
 from ..utils.extmath import fast_dot, svd_flip
 
 
-def _weighted_average(x, y, x_weight, y_weight):
-    num = x * x_weight + y * y_weight
-    denom = x_weight + y_weight
+def _mean_update(old_mean, new_mean, old_sample_count, new_sample_count):
+    """Minibatch mean update."""
+    num = old_mean * old_sample_count + new_mean * new_sample_count
+    denom = old_sample_count + new_sample_count
     return num / denom
 
 
+def _calc_sum_and_var(X):
+    """Calculate Youngs and Cramer components (T and S)."""
+    stored_sum = np.sum(X, axis=0)
+    unnormalized_variance = np.sum((X - 1. / X.shape[0] *
+                                    stored_sum) ** 2, axis=0)
+    return stored_sum, unnormalized_variance
+
+
+def _variance_update(old_sum, new_sum, old_var, new_var, old_sample_count,
+                     new_sample_count):
+    """Youngs and Cramer minibatch update."""
+    batch_sum = old_sum + new_sum
+    n = new_sample_count
+    m = old_sample_count
+    partial_var = float(m) / (n * (m + n)) * (n / float(m) * old_sum
+                                              - new_sum) ** 2
+    batch_var = old_var + new_var + partial_var
+    return batch_sum, batch_var
+
+
 class IncrementalPCA(BaseEstimator, TransformerMixin):
     """Incremental principal components analysis (IPCA).
 
@@ -124,9 +145,6 @@ def fit(self, X, y=None):
         if hasattr(self, "components_"):
             del self.components_
             del self.mean_
-            del self.explained_variance_
-            del self.explained_variance_ratio_
-            del self._explained_variance_sum
         self.samples_seen_ = 0
         X = array2d(X)
         n_samples, n_features = X.shape
@@ -178,19 +196,19 @@ def partial_fit(self, X, y=None):
             U, S, V = linalg.svd(X, full_matrices=False)
             U, V = svd_flip(U, V, u_based_decision=False)
             components = V[:n_components]
-
-            explained_variance = (S ** 2) / n_samples
-            explained_variance_sum = np.sum(explained_variance)
-            explained_variance_ratio = (explained_variance /
-                                        explained_variance_sum)
+            stored_sum, unnormalized_variance = _calc_sum_and_var(X)
+            explained_variance = S ** 2 / n_samples
+            variance_sum = np.sum(unnormalized_variance / n_samples)
+            explained_variance_ratio = explained_variance / variance_sum
         else:
             old_components = self.components_
             old_mean = self.mean_
             old_sample_count = self.samples_seen_
             new_sample_count = n_samples
             new_mean = X.mean(axis=0)
-            mean = _weighted_average(old_mean, new_mean,
-                                     old_sample_count, new_sample_count)
+            mean = _mean_update(old_mean, new_mean, old_sample_count,
+                                new_sample_count)
+
             X -= new_mean
             append_vals = np.sqrt((old_sample_count * new_sample_count) /
                                   (old_sample_count + new_sample_count))
@@ -204,22 +222,25 @@ def partial_fit(self, X, y=None):
             U, V = svd_flip(U, V, u_based_decision=False)
             components = V[:n_components]
 
-            explained_variance = (S ** 2) / (old_sample_count + new_sample_count)
-            old_sum = self._explained_variance_sum
-            new_sum = np.sum(explained_variance)
-            var_diff = np.abs(new_sum - old_sum)
-            explained_variance_sum = old_sum + (float(new_sample_count) /
-                                                old_sample_count) * var_diff
-            explained_variance_ratio = (explained_variance /
-                                        explained_variance_sum)
+            old_stored_sum = self._stored_sum
+            old_unnormalized_variance = self._unnormalized_variance
+            new_stored_sum, new_unnormalized_variance = _calc_sum_and_var(X)
+            stored_sum, unnormalized_variance = _variance_update(
+                old_stored_sum, new_stored_sum, old_unnormalized_variance,
+                new_unnormalized_variance, old_sample_count, new_sample_count)
+            explained_variance = S ** 2 / (old_sample_count + new_sample_count)
+            variance_sum = np.sum(unnormalized_variance / (old_sample_count +
+                                                           new_sample_count))
+            explained_variance_ratio = explained_variance / variance_sum
 
         self.samples_seen_ += new_sample_count
         self.components_ = components[:n_components]
         self.singular_vals_ = S[:n_components]
+        self.mean_ = mean
         self.explained_variance_ = explained_variance[:n_components]
-        self._explained_variance_sum = explained_variance_sum
         self.explained_variance_ratio_ = explained_variance_ratio[:n_components]
-        self.mean_ = mean
+        self._stored_sum = stored_sum
+        self._unnormalized_variance = unnormalized_variance
         return self
 
     def transform(self, X):
diff --git a/sklearn/decomposition/tests/test_incremental_pca.py b/sklearn/decomposition/tests/test_incremental_pca.py
@@ -7,6 +7,8 @@
 
 from sklearn import datasets
 from sklearn.decomposition import PCA, IncrementalPCA
+from sklearn.decomposition.incremental_pca import _calc_sum_and_var
+from sklearn.decomposition.incremental_pca import _variance_update
 
 iris = datasets.load_iris()
 
@@ -144,6 +146,29 @@ def test_incremental_pca_against_pca_random_data():
     assert_almost_equal(np.abs(Y_pca), np.abs(Y_ipca), 1)
 
 
+def test_update_formulas():
+    # Doggie data from http://www.mathsisfun.com/data/standard-deviation.html
+    A = np.array([[600, 470, 170, 430, 300],
+                  [600, 470, 170, 430, 300],
+                  [600, 470, 170, 430, 300],
+                  [600, 470, 170, 430, 300]]).T
+    idx = 2
+    X1 = A[:idx, :]
+    X2 = A[idx:, :]
+
+    old_sample_count = X1.shape[0]
+    new_sample_count = X2.shape[0]
+    old_stored_sum, old_unnormalized_variance = _calc_sum_and_var(X1)
+    new_stored_sum, new_unnormalized_variance = _calc_sum_and_var(X2)
+    stored_sum, unnormalized_variance = _variance_update(
+        old_stored_sum, new_stored_sum, old_unnormalized_variance,
+        new_unnormalized_variance, old_sample_count, new_sample_count)
+    sums = np.array([1970] * A.shape[1])
+    variances = np.array([108520] * A.shape[1])
+    assert_almost_equal(stored_sum, sums, 6)
+    assert_almost_equal(unnormalized_variance, variances, 6)
+
+
 if __name__ == '__main__':
     import nose
     nose.run(argv=['', __file__])
