diff --git a/bob/learn/em/data/gmm_ML.hdf5 b/bob/learn/em/data/gmm_ML.hdf5
index 74269fe3d824aa877e59609097828ee922d5dbc4..0326c186f11ad38387b3a4c4b1e3cc7a66d32b80 100644
Binary files a/bob/learn/em/data/gmm_ML.hdf5 and b/bob/learn/em/data/gmm_ML.hdf5 differ
diff --git a/bob/learn/em/mixture/gmm.py b/bob/learn/em/mixture/gmm.py
index d64b6f190cdc387fca2fa8fbc17b5f2e82cb3f5b..4370d5b01822b8ca50350bd181be41ec2ea41327 100644
--- a/bob/learn/em/mixture/gmm.py
+++ b/bob/learn/em/mixture/gmm.py
@@ -439,7 +439,7 @@ class GMMMachine(BaseEstimator):
         try:
             version_major, version_minor = hdf5.get("meta_file_version")[()].split(".")
             logger.debug(
-                f"Reading a GMMStats HDF5 file of version {version_major}.{version_minor}"
+                f"Reading a GMMMachine HDF5 file of version {version_major}.{version_minor}"
             )
         except (TypeError, RuntimeError):
             version_major, version_minor = 0, 0
@@ -460,12 +460,12 @@ class GMMMachine(BaseEstimator):
                 update_variances=hdf5["update_variances"][()],
                 update_weights=hdf5["update_weights"][()],
             )
-            gaussians = hdf5["gaussians"]
-            self.means = gaussians["means"][()]
-            self.variances = gaussians["variances"][()]
-            self.variance_thresholds = gaussians["variance_thresholds"][()]
+            gaussians_group = hdf5["gaussians"]
+            self.means = gaussians_group["means"][()]
+            self.variances = gaussians_group["variances"][()]
+            self.variance_thresholds = gaussians_group["variance_thresholds"][()]
         else:  # Legacy file version
-            logger.info("Loading a legacy HDF5 stats file.")
+            logger.info("Loading a legacy HDF5 machine file.")
             n_gaussians = int(hdf5["m_n_gaussians"][()])
             g_means = []
             g_variances = []
@@ -475,14 +475,11 @@ class GMMMachine(BaseEstimator):
                 g_means.append(gaussian_group["m_mean"][()])
                 g_variances.append(gaussian_group["m_variance"][()])
                 g_variance_thresholds.append(gaussian_group["m_variance_thresholds"][()])
-            self = cls(
-                n_gaussians=n_gaussians,
-                ubm=ubm,
-                weights=hdf5["m_weights"][()],
-            )
-            self.means = np.array(g_means)
-            self.variances = np.array(g_variances)
-            self.variance_thresholds = np.array(g_variance_thresholds)
+            weights = hdf5["m_weights"][()].reshape(n_gaussians)
+            self = cls(n_gaussians=n_gaussians, ubm=ubm, weights=weights)
+            self.means = np.array(g_means).reshape(n_gaussians,-1)
+            self.variances = np.array(g_variances).reshape(n_gaussians,-1)
+            self.variance_thresholds = np.array(g_variance_thresholds).reshape(n_gaussians,-1)
         return self
 
     def save(self, hdf5):
diff --git a/bob/learn/em/test/test_gmm.py b/bob/learn/em/test/test_gmm.py
index 40d535a9fcaa1dba3d144769ce2be4618ace0ed1..d69415d648bf21cdac42664249df42025e096157 100644
--- a/bob/learn/em/test/test_gmm.py
+++ b/bob/learn/em/test/test_gmm.py
@@ -547,7 +547,7 @@ def test_map_transformer():
 ## Tests from `test_em.py`
 
 def loadGMM():
-    gmm = GMMMachine(2)
+    gmm = GMMMachine(n_gaussians=2)
 
     gmm.weights = bob.io.base.load(datafile("gmm.init_weights.hdf5", __name__, path="../data/"))
     gmm.means = bob.io.base.load(datafile("gmm.init_means.hdf5", __name__, path="../data/"))
@@ -578,13 +578,11 @@ def test_gmm_ML_1():
     gmm.update_weights = True
     gmm = gmm.fit(ar)
 
-    #config = HDF5File(datafile("gmm_ML.hdf5", __name__), "w")
-    #gmm.save(config)
+    # Generate reference
+    # gmm.save(HDF5File(datafile("gmm_ML.hdf5", __name__, path="../data"), "w"))
 
-    gmm_ref = GMMMachine.from_hdf5(HDF5File(datafile("gmm_ML.hdf5", __name__, path="../data"), "r")) # TODO update the ref file(s)
-    gmm_ref_32bit_debug = GMMMachine.from_hdf5(HDF5File(datafile("gmm_ML_32bit_debug.hdf5", __name__, path="../data/"), "r"))
-    gmm_ref_32bit_release = GMMMachine.from_hdf5(HDF5File(datafile("gmm_ML_32bit_release.hdf5", __name__, path="../data/"), "r"))
-    assert (gmm == gmm_ref)  # or (gmm == gmm_ref_32bit_release) or (gmm == gmm_ref_32bit_debug)
+    gmm_ref = GMMMachine.from_hdf5(HDF5File(datafile("gmm_ML.hdf5", __name__, path="../data"), "r"))
+    assert gmm == gmm_ref
 
 
 def test_gmm_ML_2():
@@ -592,7 +590,7 @@ def test_gmm_ML_2():
     ar = bob.io.base.load(datafile("dataNormalized.hdf5", __name__, path="../data/"))
 
     # Initialize GMMMachine
-    gmm = GMMMachine(5, 45)
+    gmm = GMMMachine(n_gaussians=5)
     gmm.means = bob.io.base.load(datafile("meansAfterKMeans.hdf5", __name__, path="../data/")).astype("float64")
     gmm.variances = bob.io.base.load(datafile("variancesAfterKMeans.hdf5", __name__, path="../data/")).astype("float64")
     gmm.weights = np.exp(bob.io.base.load(datafile("weightsAfterKMeans.hdf5", __name__, path="../data/")).astype("float64"))
@@ -628,7 +626,7 @@ def test_gmm_ML_parallel():
     ar = da.array(bob.io.base.load(datafile("dataNormalized.hdf5", __name__, path="../data/")))
 
     # Initialize GMMMachine
-    gmm = GMMMachine(5, 45)
+    gmm = GMMMachine(n_gaussians=5)
     gmm.means = bob.io.base.load(datafile("meansAfterKMeans.hdf5", __name__, path="../data/")).astype("float64")
     gmm.variances = bob.io.base.load(datafile("variancesAfterKMeans.hdf5", __name__, path="../data/")).astype("float64")
     gmm.weights = np.exp(bob.io.base.load(datafile("weightsAfterKMeans.hdf5", __name__, path="../data/")).astype("float64"))
@@ -703,6 +701,8 @@ def test_gmm_MAP_2():
     gmm.variances = variances
     gmm.weights = weights
 
+    gmm = gmm.fit(data)
+
     gmm_adapted = GMMMachine(
         n_gaussians=2,
         trainer="map",
@@ -717,7 +717,6 @@ def test_gmm_MAP_2():
     gmm_adapted.variances = variances
     gmm_adapted.weights = weights
 
-    gmm = gmm.fit(data)
 
     gmm_adapted = gmm_adapted.fit(data)