Merge pull request scikit-learn#6201 from GaelVaroquaux/more_ledoit_wold_tests

[MRG+1] TST: More tests for Ledoit-Wolf

Author: ogrisel

sklearn/covariance/shrunk_covariance_.py

@@ -198,6 +198,9 @@ def ledoit_wolf_shrinkage(X, assume_centered=False, block_size=1000):
     if not assume_centered:
         X = X - X.mean(0)
 
+    # A non-blocked version of the computation is present in the tests
+    # in tests/test_covariance.py
+
     # number of blocks to split the covariance matrix into
     n_splits = int(n_features / block_size)
     X2 = X ** 2
@@ -232,6 +235,8 @@ def ledoit_wolf_shrinkage(X, assume_centered=False, block_size=1000):
     delta = delta_ - 2. * mu * emp_cov_trace.sum() + n_features * mu ** 2
     delta /= n_features
     # get final beta as the min between beta and delta
+    # We do this to prevent shrinking more than "1", which whould invert
+    # the value of covariances
     beta = min(beta, delta)
     # finally get shrinkage
     shrinkage = 0 if beta == 0 else beta / delta

sklearn/covariance/tests/test_covariance.py

@@ -112,6 +112,7 @@ def test_ledoit_wolf():
     lw = LedoitWolf(assume_centered=True)
     lw.fit(X_centered)
     shrinkage_ = lw.shrinkage_
+
     score_ = lw.score(X_centered)
     assert_almost_equal(ledoit_wolf_shrinkage(X_centered,
                                               assume_centered=True),
@@ -186,6 +187,38 @@ def test_ledoit_wolf():
     assert(lw.precision_ is None)
 
 
+def _naive_ledoit_wolf_shrinkage(X):
+    # A simple implementation of the formulas from Ledoit & Wolf
+
+    # The computation below achieves the following computations of the
+    # "O. Ledoit and M. Wolf, A Well-Conditioned Estimator for
+    # Large-Dimensional Covariance Matrices"
+    # beta and delta are given in the beginning of section 3.2
+    n_samples, n_features = X.shape
+    emp_cov = empirical_covariance(X, assume_centered=False)
+    mu = np.trace(emp_cov) / n_features
+    delta_ = emp_cov.copy()
+    delta_.flat[::n_features + 1] -= mu
+    delta = (delta_ ** 2).sum() / n_features
+    X2 = X ** 2
+    beta_ = 1. / (n_features * n_samples) \
+        * np.sum(np.dot(X2.T, X2) / n_samples - emp_cov ** 2)
+
+    beta = min(beta_, delta)
+    shrinkage = beta / delta
+    return shrinkage
+
+
+def test_ledoit_wolf_small():
+    # Compare our blocked implementation to the naive implementation
+    X_small = X[:, :4]
+    lw = LedoitWolf()
+    lw.fit(X_small)
+    shrinkage_ = lw.shrinkage_
+
+    assert_almost_equal(shrinkage_, _naive_ledoit_wolf_shrinkage(X_small))
+
+
 def test_ledoit_wolf_large():
     # test that ledoit_wolf doesn't error on data that is wider than block_size
     rng = np.random.RandomState(0)