Add a sample rejection procedure based on a simple likelihood threshold.

picard/_core_picard.py

@@ -8,11 +8,13 @@
 from scipy.linalg import expm
 
 from .densities import Tanh
+from ._tools import fuse_mask
 
 
 def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7,
                 max_iter=1000, tol=1e-7, lambda_min=0.01, ls_tries=10,
-                verbose=False, covariance=None):
+                verbose=False, covariance=None, ll_reject=None,
+                reject_every=None):
     '''Runs the Picard algorithm
 
     The algorithm is detailed in::
@@ -78,12 +80,25 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7,
     current_loss = _loss(Y, W, density, signs, ortho, extended)
     requested_tolerance = False
     sign_change = False
+    rejection = ll_reject is not None
     gradient_norm = 1.
     if extended:
         if covariance is None:  # Need this for extended
             covariance = X.dot(X.T) / T
         C = covariance.copy()
+    if rejection:  # create a mask
+        mask = np.zeros(T, dtype=np.bool_)
     for n in range(max_iter):
+        # Rejection
+        if rejection and (n + 1) % reject_every == 0:
+            ll = density.log_lik(Y)
+            rejected = np.sum(ll > ll_reject, axis=0) > 0
+            # Exclude rejected data
+            Y = Y[:, ~rejected]
+            # Fuse the mask to keep track of the total mask
+            mask = fuse_mask(mask, rejected)
+            # Compute the new loss
+            current_loss = _loss(Y, W, density, signs, ortho, extended)
         # Compute the score function
         psiY, psidY = density.score_and_der(Y)
         # Compute the relative gradient and the Hessian off-diagonal
@@ -169,6 +184,8 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7,
                  n_iterations=n)
     if extended:
         infos['signs'] = signs
+    if rejection:
+        infos['rejected'] = mask
     return Y, W, infos
 
 

picard/_tools.py

@@ -89,3 +89,13 @@ def _ica_par(X, fun, max_iter, w_init, verbose):
     if verbose:
         print('Running Picard...')
     return W
+
+
+def fuse_mask(mask1, mask2):
+    idx_m2 = 0
+    N = len(mask1)
+    for i in range(N):
+        if not mask1[i]:
+            mask1[i] = mask2[idx_m2]
+            idx_m2 += 1
+    return mask1

picard/solver.py

@@ -15,9 +15,9 @@
 
 def picard(X, fun='tanh', n_components=None, ortho=True, extended=None,
            whiten=True, return_X_mean=False, return_n_iter=False,
-           centering=True, max_iter=100, tol=1e-07, m=7,  ls_tries=10,
-           lambda_min=0.01, check_fun=True, w_init=None, fastica_it=None,
-           random_state=None, verbose=False):
+           centering=True, ll_reject=None, reject_every=5, max_iter=100,
+           tol=1e-07, m=7,  ls_tries=10, lambda_min=0.01, check_fun=True,
+           w_init=None, fastica_it=None, random_state=None, verbose=False):
     """Perform Independent Component Analysis.
 
     Parameters
@@ -67,6 +67,14 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None,
     centering : bool, optional
         If True, X is mean corrected.
 
+    ll_reject : None or float, optional
+        If None, no sample rejection is performed. Otherwise, samples with a
+        negative log-likelihood > reject are rejected every `reject_every`
+        iteration.
+
+    reject_every : int, optional
+        Samples are rejected every `reject_every` iteration
+
     max_iter : int, optional
         Maximum number of iterations to perform.
 
@@ -129,9 +137,15 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None,
     n_iter : int
         Number of iterations taken to converge. This is
         returned only when return_n_iter is set to `True`.
+
+    rejected : boolean array, shape (n_samples, )
+        A mask indicating which samples have been rejected (1 = rejected).
+        Returned if `ll_reject` is not None.
     """
     random_state = check_random_state(random_state)
 
+    rejection = (ll_reject is not None)
+
     # Behaves like Fastica if ortho, standard infomax if not ortho
     if extended is None:
         extended = ortho
@@ -200,7 +214,8 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None,
     kwargs = {'density': fun, 'm': m, 'max_iter': max_iter, 'tol': tol,
               'lambda_min': lambda_min, 'ls_tries': ls_tries,
               'verbose': verbose, 'ortho': ortho, 'extended': extended,
-              'covariance': covariance}
+              'covariance': covariance, 'll_reject': ll_reject,
+              'reject_every': reject_every}
 
     with np.errstate(all='raise'):  # So code breaks if overflow/Nans
         Y, W, infos = core_picard(X1, **kwargs)
@@ -217,6 +232,24 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None,
     if not whiten:
         K = None
     n_iter = infos['n_iterations']
+    if rejection:
+        rejected = infos['rejected']
+        if return_X_mean:
+            if centering:
+                if return_n_iter:
+                    return K, W, Y, X_mean, n_iter, rejected
+                else:
+                    return K, W, Y, X_mean, rejected
+            else:
+                if return_n_iter:
+                    return K, W, Y, np.zeros(p), n_iter, rejected
+                else:
+                    return K, W, Y, np.zeros(p), rejected
+        else:
+            if return_n_iter:
+                return K, W, Y, n_iter, rejected
+            else:
+                return K, W, Y, rejected
     if return_X_mean:
         if centering:
             if return_n_iter:

picard/tests/test_solver.py

@@ -96,6 +96,31 @@ def test_pre_fastica():
                         err_msg=err_msg)
 
 
+def test_rejection():
+    N, T = 3, 1000
+    rng = np.random.RandomState(42)
+    fun = 'tanh'
+    extended = False
+    A = rng.randn(N, N)
+    S = rng.laplace(size=(N, T))
+    X = np.dot(A, S)
+    # Add some strong artifacts
+    X[1, 10] = 1e2
+    X[2, 20] = 1e2
+    X[0, 20] = -1e2
+    X[2, 100] = 1e2
+    K, W, Y, rejected = picard(X, fun=fun, ortho=False, extended=extended,
+                               random_state=0, ll_reject=15)
+    for loc in [10, 20, 100]:
+        assert rejected[loc]
+    assert np.sum(rejected) == 3
+    assert_equal(W.shape, A.shape)
+    assert_equal(K.shape, A.shape)
+    WA = W.dot(K).dot(A)
+    WA = permute(WA)  # Permute and scale
+    assert_allclose(WA, np.eye(N), rtol=0, atol=1e-1)
+
+
 def test_picard():
     N, T = 3, 1000
     rng = np.random.RandomState(42)