diff --git a/bob/learn/em/factor_analysis.py b/bob/learn/em/factor_analysis.py
index 8a76d29d084e68abc9d3cc1dbe4f9bcb5eaecbb0..18cdaa0b669fe0179afb8aac80eb9ec1703a2efd 100644
--- a/bob/learn/em/factor_analysis.py
+++ b/bob/learn/em/factor_analysis.py
@@ -89,9 +89,12 @@ class FactorAnalysisBase(BaseEstimator):
         em_iterations=10,
         seed=0,
         ubm=None,
-        **gmm_kwargs,
+        gmm_kwargs=None,
+        **kwargs,
     ):
+        super().__init__(**kwargs)
         self.ubm = ubm
+        self.gmm_kwargs = gmm_kwargs
         self.em_iterations = em_iterations
         self.seed = seed
 
@@ -1146,7 +1149,7 @@ class FactorAnalysisBase(BaseEstimator):
 
         """
 
-        return self.score(model, self.ubm.acc_statistics(data))
+        return self.score(model, self.ubm.transform(data))
 
 
 class ISVMachine(FactorAnalysisBase):
@@ -1338,7 +1341,7 @@ class ISVMachine(FactorAnalysisBase):
             z
 
         """
-        return self.enroll([self.ubm.acc_statistics(X)], iterations)
+        return self.enroll([self.ubm.transform(X)], iterations)
 
     def score(self, latent_z, data):
         """
@@ -1764,7 +1767,7 @@ class JFAMachine(FactorAnalysisBase):
             z
 
         """
-        return self.enroll([self.ubm.acc_statistics(X)], iterations)
+        return self.enroll([self.ubm.transform(X)], iterations)
 
     def fit(self, X, y):
         """
diff --git a/bob/learn/em/test/test_jfa.py b/bob/learn/em/test/test_jfa.py
index de4f76115a8611bce13955d719c082a35eb4cc28..8e20b4a4f965cbb31f59caa01c36e876c76752f0 100644
--- a/bob/learn/em/test/test_jfa.py
+++ b/bob/learn/em/test/test_jfa.py
@@ -16,35 +16,26 @@ def test_JFAMachine():
     eps = 1e-10
 
     # Creates a UBM
-    weights = np.array([0.4, 0.6], "float64")
-    means = np.array([[1, 6, 2], [4, 3, 2]], "float64")
-    variances = np.array([[1, 2, 1], [2, 1, 2]], "float64")
     ubm = GMMMachine(2, 3)
-    ubm.weights = weights
-    ubm.means = means
-    ubm.variances = variances
+    ubm.weights = np.array([0.4, 0.6], "float64")
+    ubm.means = np.array([[1, 6, 2], [4, 3, 2]], "float64")
+    ubm.variances = np.array([[1, 2, 1], [2, 1, 2]], "float64")
 
     # Defines GMMStats
     gs = GMMStats(2, 3)
-    log_likelihood = -3.0
-    T = 1
-    n = np.array([0.4, 0.6], "float64")
-    sumpx = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], "float64")
-    sumpxx = np.array([[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]], "float64")
-    gs.log_likelihood = log_likelihood
-    gs.t = T
-    gs.n = n
-    gs.sum_px = sumpx
-    gs.sum_pxx = sumpxx
+    gs.log_likelihood = -3.0
+    gs.t = 1
+    gs.n = np.array([0.4, 0.6], "float64")
+    gs.sum_px = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], "float64")
+    gs.sum_pxx = np.array([[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]], "float64")
 
     # Creates a JFAMachine
-    U = np.array([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], "float64")
-    V = np.array([[6, 5], [4, 3], [2, 1], [1, 2], [3, 4], [5, 6]], "float64")
-    d = np.array([0, 1, 0, 1, 0, 1], "float64")
     m = JFAMachine(ubm, 2, 2, em_iterations=10)
-    m.U = U
-    m.V = V
-    m.D = d
+    m.U = np.array(
+        [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], "float64"
+    )
+    m.V = np.array([[6, 5], [4, 3], [2, 1], [1, 2], [3, 4], [5, 6]], "float64")
+    m.D = np.array([0, 1, 0, 1, 0, 1], "float64")
 
     # Preparing the model
     y = np.array([1, 2], "float64")
@@ -53,15 +44,14 @@ def test_JFAMachine():
 
     score_ref = -2.111577181208289
     score = m.score(model, gs)
-    assert abs(score_ref - score) < eps
+    np.testing.assert_allclose(score, score_ref, atol=eps)
 
     # Scoring with numpy array
     np.random.seed(0)
     X = np.random.normal(loc=0.0, scale=1.0, size=(50, 3))
     score_ref = 2.028009315286946
     score = m.score_with_array(model, X)
-
-    assert abs(score_ref - score) < eps
+    np.testing.assert_allclose(score, score_ref, atol=eps)
 
 
 def test_ISVMachine():
@@ -69,48 +59,36 @@ def test_ISVMachine():
     eps = 1e-10
 
     # Creates a UBM
-    weights = np.array([0.4, 0.6], "float64")
-    means = np.array([[1, 6, 2], [4, 3, 2]], "float64")
-    variances = np.array([[1, 2, 1], [2, 1, 2]], "float64")
     ubm = GMMMachine(2, 3)
-    ubm.weights = weights
-    ubm.means = means
-    ubm.variances = variances
+    ubm.weights = np.array([0.4, 0.6], "float64")
+    ubm.means = np.array([[1, 6, 2], [4, 3, 2]], "float64")
+    ubm.variances = np.array([[1, 2, 1], [2, 1, 2]], "float64")
 
     # Creates a ISVMachine
-    U = np.array([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], "float64")
-    # V = numpy.array([[0], [0], [0], [0], [0], [0]], 'float64')
-    d = np.array([0, 1, 0, 1, 0, 1], "float64")
     isv_machine = ISVMachine(ubm, r_U=2, em_iterations=10)
-    isv_machine.U = U
-    # base.v = V
-    isv_machine.D = d
+    isv_machine.U = np.array(
+        [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], "float64"
+    )
+    # base.v = numpy.array([[0], [0], [0], [0], [0], [0]], 'float64')
+    isv_machine.D = np.array([0, 1, 0, 1, 0, 1], "float64")
 
     # Defines GMMStats
     gs = GMMStats(2, 3)
-    log_likelihood = -3.0
-    T = 1
-    n = np.array([0.4, 0.6], "float64")
-    sumpx = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], "float64")
-    sumpxx = np.array([[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]], "float64")
-    gs.log_likelihood = log_likelihood
-    gs.t = T
-    gs.n = n
-    gs.sum_px = sumpx
-    gs.sum_pxx = sumpxx
+    gs.log_likelihood = -3.0
+    gs.t = 1
+    gs.n = np.array([0.4, 0.6], "float64")
+    gs.sum_px = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], "float64")
+    gs.sum_pxx = np.array([[10.0, 20.0, 30.0], [40.0, 50.0, 60.0]], "float64")
 
     # Enrolled model
     latent_z = np.array([3, 4, 1, 2, 0, 1], "float64")
     score = isv_machine.score(latent_z, gs)
     score_ref = -3.280498193082100
-
-    assert abs(score_ref - score) < eps
+    np.testing.assert_allclose(score, score_ref, atol=eps)
 
     # Scoring with numpy array
     np.random.seed(0)
     X = np.random.normal(loc=0.0, scale=1.0, size=(50, 3))
     score_ref = -1.2343813195374242
-
     score = isv_machine.score_with_array(latent_z, X)
-
-    assert abs(score_ref - score) < eps
+    np.testing.assert_allclose(score, score_ref, atol=eps)
diff --git a/bob/learn/em/test/test_jfa_trainer.py b/bob/learn/em/test/test_jfa_trainer.py
index f2b08d1d9ba055fb918298c783059b3e5c43c72f..cc3d6ad1b57c22fe847f3f26e2674831f49b5c0b 100644
--- a/bob/learn/em/test/test_jfa_trainer.py
+++ b/bob/learn/em/test/test_jfa_trainer.py
@@ -163,9 +163,9 @@ def test_JFATrainAndEnrol():
     )
 
     eps = 1e-10
-    assert np.allclose(it.V, v_ref, eps)
-    assert np.allclose(it.U, u_ref, eps)
-    assert np.allclose(it.D, d_ref, eps)
+    np.testing.assert_allclose(it.V, v_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.U, u_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.D, d_ref, rtol=eps, atol=1e-8)
 
     # Calls the enroll function
 
@@ -209,8 +209,8 @@ def test_JFATrainAndEnrol():
         "float64",
     )
 
-    assert np.allclose(latent_y, y_ref, eps)
-    assert np.allclose(latent_z, z_ref, eps)
+    np.testing.assert_allclose(latent_y, y_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(latent_z, z_ref, rtol=eps, atol=1e-8)
 
 
 def test_JFATrainAndEnrolWithNumpy():
@@ -262,9 +262,9 @@ def test_JFATrainAndEnrolWithNumpy():
     )
 
     eps = 1e-10
-    assert np.allclose(it.V, v_ref, eps)
-    assert np.allclose(it.U, u_ref, eps)
-    assert np.allclose(it.D, d_ref, eps)
+    np.testing.assert_allclose(it.V, v_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.U, u_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.D, d_ref, rtol=eps, atol=1e-8)
 
     """
     Calls the enroll function with arrays as input
@@ -275,20 +275,18 @@ def test_JFATrainAndEnrolWithNumpy():
     latent_y_ref = np.array([-0.13922039, 0.10686916])
     latent_z_ref = np.array(
         [
-            [
-                -1.37073043e-17,
-                1.15641870e-12,
-                -8.29922598e-10,
-                -4.17108194e-16,
-                -2.27107305e-09,
-                2.94293314e-18,
-            ]
+            -1.37073043e-17,
+            1.15641870e-12,
+            -8.29922598e-10,
+            -4.17108194e-16,
+            -2.27107305e-09,
+            2.94293314e-18,
         ]
     )
 
     latent_y, latent_z = it.enroll_with_array(X)
-    assert np.allclose(latent_z, latent_z_ref, eps)
-    assert np.allclose(latent_y, latent_y_ref, eps)
+    np.testing.assert_allclose(latent_z, latent_z_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(latent_y, latent_y_ref, rtol=eps, atol=1e-8)
 
 
 def test_ISVTrainAndEnrol():
@@ -319,12 +317,14 @@ def test_ISVTrainAndEnrol():
     )
     z_ref = np.array(
         [
-            -0.079315777443826,
-            0.092702428248543,
-            -0.342488761656616,
-            -0.059922635809136,
-            0.133539981073604,
-            0.213118695516570,
+            [
+                -0.079315777443826,
+                0.092702428248543,
+                -0.342488761656616,
+                -0.059922635809136,
+                0.133539981073604,
+                0.213118695516570,
+            ]
         ],
         "float64",
     )
@@ -346,8 +346,8 @@ def test_ISVTrainAndEnrol():
     it.U = copy.deepcopy(M_u)
     it = it.fit(TRAINING_STATS_X, TRAINING_STATS_y)
 
-    assert np.allclose(it.D, d_ref, eps)
-    assert np.allclose(it.U, u_ref, eps)
+    np.testing.assert_allclose(it.D, d_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.U, u_ref, rtol=eps, atol=1e-8)
 
     """
     Calls the enroll function
@@ -379,7 +379,7 @@ def test_ISVTrainAndEnrol():
 
     gse = [gse1, gse2]
     latent_z = it.enroll(gse, 5)
-    assert np.allclose(latent_z, z_ref, eps)
+    np.testing.assert_allclose(latent_z, z_ref, rtol=eps, atol=1e-8)
 
 
 def test_ISVTrainAndEnrolWithNumpy():
@@ -426,8 +426,8 @@ def test_ISVTrainAndEnrolWithNumpy():
     it.U = copy.deepcopy(M_u)
     it = it.fit(TRAINING_STATS_X, TRAINING_STATS_y)
 
-    assert np.allclose(it.D, d_ref, eps)
-    assert np.allclose(it.U, u_ref, eps)
+    np.testing.assert_allclose(it.D, d_ref, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(it.U, u_ref, rtol=eps, atol=1e-8)
 
     """
     Calls the enroll function with arrays as input
@@ -449,7 +449,7 @@ def test_ISVTrainAndEnrolWithNumpy():
     )
 
     latent_z = it.enroll_with_array(X)
-    assert np.allclose(latent_z, latent_z_ref, eps)
+    np.testing.assert_allclose(latent_z, latent_z_ref, rtol=eps, atol=1e-8)
 
 
 def test_JFATrainInitialize():
@@ -477,9 +477,9 @@ def test_JFATrainInitialize():
     v2 = it.V
     d2 = it.D
 
-    assert np.allclose(u1, u2, eps)
-    assert np.allclose(v1, v2, eps)
-    assert np.allclose(d1, d2, eps)
+    np.testing.assert_allclose(u1, u2, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(v1, v2, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(d1, d2, rtol=eps, atol=1e-8)
 
 
 def test_ISVTrainInitialize():
@@ -505,5 +505,5 @@ def test_ISVTrainInitialize():
     u2 = it.U
     d2 = it.D
 
-    assert np.allclose(u1, u2, eps)
-    assert np.allclose(d1, d2, eps)
+    np.testing.assert_allclose(u1, u2, rtol=eps, atol=1e-8)
+    np.testing.assert_allclose(d1, d2, rtol=eps, atol=1e-8)