diff --git a/bob/bio/gmm/algorithm/GMM.py b/bob/bio/gmm/algorithm/GMM.py
index 672ed2309b386db0571085f23f3b2aca7ea3d230..3d9d23eadbe0bd8c4902db8b973a880ecc44a692 100644
--- a/bob/bio/gmm/algorithm/GMM.py
+++ b/bob/bio/gmm/algorithm/GMM.py
@@ -230,9 +230,8 @@ class GMM(BioAlgorithm, BaseEstimator):
 
         logger.debug(f"scoring {biometric_reference}, {probe}")
         if not isinstance(probe, GMMStats):
-            probe = self.project(
-                probe
-            )  # Projection is done here instead of transform (or it would be applied to enrollment data too...)
+            # Projection is done here instead of in transform (or it would be applied to enrollment data too...)
+            probe = self.project(probe)
         return self.scoring_function(
             models_means=[biometric_reference],
             ubm=self.ubm,
@@ -265,26 +264,22 @@ class GMM(BioAlgorithm, BaseEstimator):
             ubm=self.ubm,
             test_stats=stats,
             frame_length_normalization=True,
-        )
+        ).reshape((-1,))
 
-    def score_for_multiple_probes(self, model, probes):
+    def score_for_multiple_probes(self, biometric_reference, probes):
         """This function computes the score between the given model and several given probe files."""
-        logger.debug(f"scoring {model}, {probes}")
-        assert isinstance(model, GMMMachine)
+        logger.debug(f"scoring {biometric_reference}, {probes}")
+        assert isinstance(biometric_reference, GMMMachine)
         stats = [
             self.project(probe) if not isinstance(probe, GMMStats) else probe
             for probe in probes
         ]
-        return (
-            self.scoring_function(
-                models_means=model.means,
-                ubm=self.ubm,
-                test_stats=stats,
-                frame_length_normalization=True,
-            )
-            .mean()
-            .reshape((-1,))
-        )
+        return self.scoring_function(
+            models_means=biometric_reference.means,
+            ubm=self.ubm,
+            test_stats=stats,
+            frame_length_normalization=True,
+        ).reshape((-1,))
 
     def fit(self, X, y=None, **kwargs):
         """Trains the UBM."""
diff --git a/bob/bio/gmm/test/data/gmm_enrolled.hdf5 b/bob/bio/gmm/test/data/gmm_enrolled.hdf5
index 2e6e337f592e53f56806bce07dd5556c3176e6ae..6466ed0a3fc524a6ab698f09540aa908f5108c22 100644
Binary files a/bob/bio/gmm/test/data/gmm_enrolled.hdf5 and b/bob/bio/gmm/test/data/gmm_enrolled.hdf5 differ
diff --git a/bob/bio/gmm/test/data/gmm_projected.hdf5 b/bob/bio/gmm/test/data/gmm_projected.hdf5
index 84437324be483253b54aff3fc2d0f1c2e1e620de..fc5e0a7c0b8f41d3b7d9a03d07c8e277f8bb5386 100644
Binary files a/bob/bio/gmm/test/data/gmm_projected.hdf5 and b/bob/bio/gmm/test/data/gmm_projected.hdf5 differ
diff --git a/bob/bio/gmm/test/data/gmm_ubm.hdf5 b/bob/bio/gmm/test/data/gmm_ubm.hdf5
index 99df686343ee83dc8332c56561675cb91fd98ad8..6a8abe718f958351135e5992621a733016127ccc 100644
Binary files a/bob/bio/gmm/test/data/gmm_ubm.hdf5 and b/bob/bio/gmm/test/data/gmm_ubm.hdf5 differ
diff --git a/bob/bio/gmm/test/test_gmm.py b/bob/bio/gmm/test/test_gmm.py
index 404becce825f821e4637c44bb5effd8cdb314bd4..e2434f0893f006453f65046a8999fa7d95ebabba 100644
--- a/bob/bio/gmm/test/test_gmm.py
+++ b/bob/bio/gmm/test/test_gmm.py
@@ -18,8 +18,10 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 import logging
+import os
 import tempfile
 
+import numpy
 import pkg_resources
 
 import bob.bio.gmm
@@ -46,6 +48,8 @@ def test_class():
         gmm1, bob.bio.base.pipelines.vanilla_biometrics.abstract_classes.BioAlgorithm
     )
     assert gmm1.number_of_gaussians == 512
+    assert "bob_fit_supports_dask_array" in gmm1._get_tags()
+    assert gmm1.transform(None) is None
 
 
 def test_training():
@@ -125,11 +129,17 @@ def test_enroll():
         "bob.bio.gmm.test", "data/gmm_enrolled.hdf5"
     )
     if regenerate_refs:
-        biometric_reference.save(reference_file)
+        gmm1.write_biometric_reference(biometric_reference, reference_file)
 
+    # Compare to pre-generated file
     gmm2 = gmm1.read_biometric_reference(reference_file)
     assert biometric_reference.is_similar_to(gmm2)
 
+    with tempfile.NamedTemporaryFile(prefix="bob_", suffix="_bioref.hdf5") as fd:
+        temp_file = fd.name
+        gmm1.write_biometric_reference(biometric_reference, reference_file)
+        assert os.path.exists(temp_file)
+
 
 def test_score():
     gmm1 = GMM(number_of_gaussians=2)
@@ -143,18 +153,27 @@ def test_score():
     probe = GMMStats.from_hdf5(
         pkg_resources.resource_filename("bob.bio.gmm.test", "data/gmm_projected.hdf5")
     )
+    probe_data = utils.random_array((20, 45), -5.0, 5.0, seed=84)
+
+    reference_score = -0.098980
 
-    reference_score = 0.045073
-    assert (
-        abs(gmm1.score(biometric_reference, probe) - reference_score) < 1e-5
-    ), "The scores differ: %3.8f, %3.8f" % (
-        gmm1.score(biometric_reference, probe),
-        reference_score,
+    numpy.testing.assert_almost_equal(
+        gmm1.score(biometric_reference, probe), reference_score, decimal=5
     )
-    assert (
-        abs(
-            gmm1.score_for_multiple_probes(biometric_reference, [probe, probe])
-            - reference_score
-        )
-        < 1e-5
+
+    multi_probes = gmm1.score_for_multiple_probes(
+        biometric_reference, [probe, probe, probe]
+    )
+    assert multi_probes.shape == (3,), multi_probes.shape
+    numpy.testing.assert_almost_equal(multi_probes, reference_score, decimal=5)
+
+    multi_refs = gmm1.score_multiple_biometric_references(
+        [biometric_reference, biometric_reference, biometric_reference], probe
+    )
+    assert multi_refs.shape == (3,), multi_refs.shape
+    numpy.testing.assert_almost_equal(multi_refs, reference_score, decimal=5)
+
+    # With not projected data
+    numpy.testing.assert_almost_equal(
+        gmm1.score(biometric_reference, probe_data), reference_score, decimal=5
     )