diff --git a/skpro/distributions/base.py b/skpro/distributions/base.py
index f6c35677..1523b17c 100644
--- a/skpro/distributions/base.py
+++ b/skpro/distributions/base.py
@@ -252,7 +252,7 @@ def cdf(self, x):
         spl = self.sample(N)
         ind = splx <= spl
 
-        return ind.groupby(level=1).mean()
+        return ind.groupby(level=1, sort=False).mean()
 
     def ppf(self, p):
         """Quantile function = percent point function = inverse cdf."""
@@ -332,7 +332,8 @@ def energy(self, x=None):
 
         # approx E[abs(X-Y)] via mean of samples of abs(X-Y) obtained from splx, sply
         spl = splx - sply
-        energy = spl.apply(np.linalg.norm, axis=1, ord=1).groupby(level=1).mean()
+        energy = spl.apply(np.linalg.norm, axis=1, ord=1)
+        energy = energy.groupby(level=1, sort=False).mean()
         energy = pd.DataFrame(energy, index=self.index, columns=["energy"])
         return energy
 
@@ -355,7 +356,7 @@ def mean(self):
         warn(self._method_error_msg("mean", fill_in=approx_method))
 
         spl = self.sample(approx_spl_size)
-        return spl.groupby(level=1).mean()
+        return spl.groupby(level=1, sort=False).mean()
 
     def var(self):
         r"""Return element/entry-wise variance of the distribution.
@@ -378,7 +379,7 @@ def var(self):
         spl1 = self.sample(approx_spl_size)
         spl2 = self.sample(approx_spl_size)
         spl = (spl1 - spl2) ** 2
-        return spl.groupby(level=1).mean()
+        return spl.groupby(level=1, sort=False).mean()
 
     def pdfnorm(self, a=2):
         r"""a-norm of pdf, defaults to 2-norm.
@@ -410,7 +411,7 @@ def pdfnorm(self, a=2):
 
         # uses formula int p(x)^a dx = E[p(X)^{a-1}], and MC approximates the RHS
         spl = [self.pdf(self.sample()) ** (a - 1) for _ in range(approx_spl_size)]
-        return pd.concat(spl, axis=0).groupby(level=1).mean()
+        return pd.concat(spl, axis=0).groupby(level=1, sort=False).mean()
 
     def _coerce_to_self_index_df(self, x):
         x = np.array(x)
diff --git a/skpro/distributions/empirical.py b/skpro/distributions/empirical.py
index 45718541..7959c808 100644
--- a/skpro/distributions/empirical.py
+++ b/skpro/distributions/empirical.py
@@ -197,9 +197,9 @@ def mean(self):
         """
         spl = self.spl
         if self.weights is None:
-            mean_df = spl.groupby(level=-1).mean()
+            mean_df = spl.groupby(level=-1, sort=False).mean()
         else:
-            mean_df = spl.groupby(level=-1).apply(
+            mean_df = spl.groupby(level=-1, sort=False).apply(
                 lambda x: np.average(x, weights=self.weights.loc[x.index], axis=0)
             )
             mean_df = pd.DataFrame(mean_df.tolist(), index=mean_df.index)
@@ -221,11 +221,11 @@ def var(self):
         spl = self.spl
         N = self._N
         if self.weights is None:
-            var_df = spl.groupby(level=-1).var(ddof=0)
+            var_df = spl.groupby(level=-1, sort=False).var(ddof=0)
         else:
             mean = self.mean()
             means = pd.concat([mean] * N, axis=0, keys=self._spl_instances)
-            var_df = spl.groupby(level=-1).apply(
+            var_df = spl.groupby(level=-1, sort=False).apply(
                 lambda x: np.average(
                     (x - means.loc[x.index]) ** 2,
                     weights=self.weights.loc[x.index],