scikit-learn · amueller · Jun 21, 2017 · Jun 11, 2017 · Jun 11, 2017 · Jun 11, 2017
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
@@ -372,6 +372,12 @@ Bug fixes
    - Add ``shuffle`` parameter to :func:`model_selection.train_test_split`.
      :issue:`#8845` by  :user:`themrmax <themrmax>`
 
+   - Fixed the implementation of `explained_variance_`
+     in :class:`decomposition.PCA`,
+     :class:`decomposition.RandomizedPCA` and
+     :class:`decomposition.IncrementalPCA`.
+     :issue:`9105` by `Hanmin Qin <https://github.com/qinhanmin2014>`_. 
+
 API changes summary
 -------------------
 

diff --git a/sklearn/decomposition/incremental_pca.py b/sklearn/decomposition/incremental_pca.py
@@ -251,7 +251,7 @@ def partial_fit(self, X, y=None, check_input=True):
 
         U, S, V = linalg.svd(X, full_matrices=False)
         U, V = svd_flip(U, V, u_based_decision=False)
-        explained_variance = S ** 2 / n_total_samples
+        explained_variance = S ** 2 / (n_total_samples - 1)
         explained_variance_ratio = S ** 2 / np.sum(col_var * n_total_samples)
 
         self.n_samples_seen_ = n_total_samples

diff --git a/sklearn/decomposition/pca.py b/sklearn/decomposition/pca.py
@@ -285,12 +285,6 @@ class PCA(_BasePCA):
     >>> print(pca.singular_values_)  # doctest: +ELLIPSIS
     [ 6.30061...]
 
-    Notes
-    -----
-    PCA uses the maximum likelihood estimate of the eigenvalues, which does not
-    include the Bessel correction, though in practice this should rarely make a
-    difference in a machine learning context.
-
     See also
     --------
     KernelPCA
@@ -346,7 +340,7 @@ def fit_transform(self, X, y=None):
 
         if self.whiten:
             # X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples)
-            U *= sqrt(X.shape[0])
+            U *= sqrt(X.shape[0] - 1)
         else:
             # X_new = X * V = U * S * V^T * V = U * S
             U *= S[:self.n_components_]
@@ -416,7 +410,7 @@ def _fit_full(self, X, n_components):
         components_ = V
 
         # Get variance explained by singular values
-        explained_variance_ = (S ** 2) / n_samples
+        explained_variance_ = (S ** 2) / (n_samples - 1)
         total_var = explained_variance_.sum()
         explained_variance_ratio_ = explained_variance_ / total_var
         singular_values_ = S.copy()  # Store the singular values.
@@ -495,8 +489,8 @@ def _fit_truncated(self, X, n_components, svd_solver):
         self.n_components_ = n_components
 
         # Get variance explained by singular values
-        self.explained_variance_ = (S ** 2) / n_samples
-        total_var = np.var(X, axis=0)
+        self.explained_variance_ = (S ** 2) / (n_samples - 1)
+        total_var = np.var(X, ddof=1, axis=0)
         self.explained_variance_ratio_ = \
             self.explained_variance_ / total_var.sum()
         self.singular_values_ = S.copy()  # Store the singular values.
@@ -714,8 +708,8 @@ def _fit(self, X):
                                  n_iter=self.iterated_power,
                                  random_state=random_state)
 
-        self.explained_variance_ = exp_var = (S ** 2) / n_samples
-        full_var = np.var(X, axis=0).sum()
+        self.explained_variance_ = exp_var = (S ** 2) / (n_samples - 1)
+        full_var = np.var(X, ddof=1, axis=0).sum()
         self.explained_variance_ratio_ = exp_var / full_var
         self.singular_values_ = S  # Store the singular values.
 

diff --git a/sklearn/decomposition/tests/test_pca.py b/sklearn/decomposition/tests/test_pca.py
@@ -174,7 +174,8 @@ def test_whitening():
         X_whitened2 = pca.transform(X_)
         assert_array_almost_equal(X_whitened, X_whitened2)
 
-        assert_almost_equal(X_whitened.std(axis=0), np.ones(n_components),
+        assert_almost_equal(X_whitened.std(ddof=1, axis=0),
+                            np.ones(n_components),
                             decimal=6)
         assert_almost_equal(X_whitened.mean(axis=0), np.zeros(n_components))
 
@@ -213,17 +214,25 @@ def test_explained_variance():
                               rpca.explained_variance_ratio_, 1)
 
     # compare to empirical variances
+    expected_result = np.linalg.eig(np.cov(X, rowvar=False))[0]
+    expected_result = sorted(expected_result, reverse=True)[:2]
+
     X_pca = pca.transform(X)
     assert_array_almost_equal(pca.explained_variance_,
-                              np.var(X_pca, axis=0))
+                              np.var(X_pca, ddof=1, axis=0))
+    assert_array_almost_equal(pca.explained_variance_, expected_result)
 
     X_pca = apca.transform(X)
     assert_array_almost_equal(apca.explained_variance_,
-                              np.var(X_pca, axis=0))
+                              np.var(X_pca, ddof=1, axis=0))
+    assert_array_almost_equal(apca.explained_variance_, expected_result)
 
     X_rpca = rpca.transform(X)
-    assert_array_almost_equal(rpca.explained_variance_, np.var(X_rpca, axis=0),
+    assert_array_almost_equal(rpca.explained_variance_,
+                              np.var(X_rpca, ddof=1, axis=0),
                               decimal=1)
+    assert_array_almost_equal(rpca.explained_variance_,
+                              expected_result, decimal=1)
 
     # Same with correlated data
     X = datasets.make_classification(n_samples, n_features,