Repository: sktime/sktime
Pull request #5714 (2024-01-09 12:00:47)

[ENH] refactor tests with parallelization backend fixtures to programmatic backend fixture lookup
This PR refactors tests with parallelization backend fixtures to programmatic backend fixture lookup.

Where previously separate locations of backend fixtures were curated, they are now all pointed towards a single interface point, `sktime.utils.parallel._get_parallel_test_fixtures`.

fkiraly: FYI @yarnabrina - unclear to me why this is failing:
![image](https://github.com/sktime/sktime/assets/7985502/5f844590-8542-4125-933c-473f3f361ee0)

fkiraly: seems to be GitHub related downtime

Patch:
diff --git a/docs/source/api_reference/alignment.rst b/docs/source/api_reference/alignment.rst
index 5306ed44c..eecec7e7e 100644
--- a/docs/source/api_reference/alignment.rst
+++ b/docs/source/api_reference/alignment.rst
@@ -6,6 +6,12 @@ Time series alignment
 The :mod:`sktime.alignment` module contains time series aligners, such as
 dynamic time warping aligners.
 
+All time series aligners in ``sktime`` can be listed using the
+``sktime.registry.all_estimators`` utility,
+using ``estimator_types="aligner"``, optionally filtered by tags.
+Valid tags can be listed using ``sktime.registry.all_tags``.
+
+
 Naive aligners
 --------------
 
diff --git a/docs/source/api_reference/param_est.rst b/docs/source/api_reference/param_est.rst
index d8d22f427..9bb1a2cdd 100644
--- a/docs/source/api_reference/param_est.rst
+++ b/docs/source/api_reference/param_est.rst
@@ -14,10 +14,6 @@ All parameter estimators in ``sktime`` can be listed using the
 using ``estimator_types="param_est"``, optionally filtered by tags.
 Valid tags can be listed using ``sktime.registry.all_tags``.
 
-.. automodule:: sktime.param_est
-    :no-members:
-    :no-inherited-members:
-
 Parameter estimators
 --------------------
 
diff --git a/docs/source/api_reference/split.rst b/docs/source/api_reference/split.rst
index 7e8228f22..493ae07c3 100644
--- a/docs/source/api_reference/split.rst
+++ b/docs/source/api_reference/split.rst
@@ -21,6 +21,8 @@ Forecasting users interested in performance evaluation are advised
 to use full backtesting instead of a single split, e.g., via ``evaluate``,
 see :ref:`forecasting API reference <forecasting_ref>`.
 
+.. currentmodule:: sktime.split
+
 .. autosummary::
     :toctree: auto_generated/
     :template: function.rst
diff --git a/pyproject.toml b/pyproject.toml
index 17aa0bdca..6e45ddd20 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -200,7 +200,7 @@ regression = [
 transformations = [
   'esig<0.10,>=0.9.7; python_version < "3.11"',
   "filterpy<1.5,>=1.4.5",
-  "holidays>=0.29,<0.40",
+  "holidays>=0.29,<0.41",
   "mne<1.6,>=1.5",
   'numba<0.59,>=0.53; python_version < "3.12"',
   "pycatch22>=0.4,<0.4.4",
diff --git a/sktime/base/_base.py b/sktime/base/_base.py
index 677b35257..f9d4c9c66 100644
--- a/sktime/base/_base.py
+++ b/sktime/base/_base.py
@@ -622,12 +622,12 @@ class BaseEstimator(BaseObject):
 
         # default retrieves all self attributes ending in "_"
         # and returns them with keys that have the "_" removed
-        fitted_params = [attr for attr in dir(obj) if attr.endswith("_")]
-        fitted_params = [x for x in fitted_params if not x.startswith("_")]
-        fitted_params = [x for x in fitted_params if hasattr(obj, x)]
-        fitted_param_dict = {p[:-1]: getattr(obj, p) for p in fitted_params}
-
-        return fitted_param_dict
+        fitted_params = {
+            attr[:-1]: getattr(obj, attr)
+            for attr in dir(obj)
+            if attr.endswith("_") and not attr.startswith("_") and hasattr(obj, attr)
+        }
+        return fitted_params
 
     def _get_fitted_params(self):
         """Get fitted parameters.
diff --git a/sktime/datatypes/_vectorize.py b/sktime/datatypes/_vectorize.py
index f8755d435..fb1a96154 100644
--- a/sktime/datatypes/_vectorize.py
+++ b/sktime/datatypes/_vectorize.py
@@ -61,9 +61,16 @@ class VectorizedDF:
     SERIES_SCITYPES = ["Series", "Panel", "Hierarchical"]
 
     def __init__(
-        self, X, y=None, iterate_as="Series", is_scitype="Panel", iterate_cols=False
+        self,
+        X,
+        y=None,
+        iterate_as="Series",
+        is_scitype="Panel",
+        iterate_cols=False,
+        remember_data=True,
     ):
-        self.X = X
+        if remember_data:
+            self.X = X
 
         if is_scitype is None:
             _, _, metadata = check_is_scitype(
@@ -89,9 +96,15 @@ class VectorizedDF:
         self._check_iterate_cols(iterate_cols)
         self.iterate_cols = iterate_cols
 
+        self.remember_data = remember_data
+
         self.converter_store = dict()
 
-        self.X_multiindex = self._init_conversion(X)
+        X_multiindex = self._init_conversion(X)
+        self.X_mi_columns = X_multiindex.columns
+        self.X_mi_index = X_multiindex.index
+        if remember_data:
+            self.X_multiindex = X_multiindex
         self.iter_indices = self._init_iter_indices()
 
         self.shape = self._iter_shape()
@@ -148,14 +161,14 @@ class VectorizedDF:
         iterate_as = self.iterate_as
         is_scitype = self.is_scitype
         iterate_cols = self.iterate_cols
-        X = self.X_multiindex
+        X_ix = self.X_mi_index
 
         if iterate_as == is_scitype:
             row_ix = None
         elif iterate_as == "Series":
-            row_ix = X.index.droplevel(-1).unique()
+            row_ix = X_ix.droplevel(-1).unique()
         elif iterate_as == "Panel":
-            row_ix = X.index.droplevel([-1, -2]).unique()
+            row_ix = X_ix.droplevel([-1, -2]).unique()
         else:
             raise RuntimeError(
                 f"unexpected value found for attribute self.iterate_as: {iterate_as}"
@@ -163,7 +176,7 @@ class VectorizedDF:
             )
 
         if iterate_cols:
-            col_ix = X.columns
+            col_ix = self.X_mi_columns
         else:
             col_ix = None
 
@@ -172,7 +185,7 @@ class VectorizedDF:
     @property
     def index(self):
         """Defaults to pandas index of X converted to pandas type."""
-        return self.X_multiindex.index
+        return self.X_mi_index
 
     def get_iter_indices(self):
         """Get indices that are iterated over in vectorization.
@@ -257,7 +270,7 @@ class VectorizedDF:
                 yield group_name, None, _enforce_index_freq(inst)
 
         iter_levels = self._iter_levels(iterate_as)
-        is_self_iter = len(iter_levels) == self.X_multiindex.index.nlevels
+        is_self_iter = len(iter_levels) == self.X_mi_index.nlevels
 
         if is_self_iter:
             yield from _iter_cols(self.X_multiindex)
@@ -285,7 +298,7 @@ class VectorizedDF:
                 iter_levels = 2
             elif iterate_as == "Series":
                 iter_levels = 1
-        return list(range(self.X_multiindex.index.nlevels - iter_levels))
+        return list(range(self.X_mi_index.nlevels - iter_levels))
 
     def _iter_shape(self, iterate_as=None, iterate_cols=None):
         """Get the number of groups and columns to iterate over.
@@ -306,11 +319,11 @@ class VectorizedDF:
             iterate_cols = self.iterate_cols
 
         iter_levels = self._iter_levels(iterate_as)
-        is_self_iter = len(iter_levels) == self.X_multiindex.index.nlevels
+        is_self_iter = len(iter_levels) == self.X_mi_index.nlevels
 
         return (
             1 if is_self_iter else self.X_multiindex.groupby(level=iter_levels).ngroups,
-            len(self.X_multiindex.columns) if iterate_cols else 1,
+            len(self.X_mi_columns) if iterate_cols else 1,
         )
 
     def as_list(self):
@@ -383,7 +396,7 @@ class VectorizedDF:
         row_ix, col_ix = self.get_iter_indices()
         force_flat = False
         if row_ix is None and col_ix is None:
-            X_mi_reconstructed = self.X_multiindex
+            X_mi_reconstructed = pd.DataFrame(df_list[0])
         elif col_ix is None:
             X_mi_reconstructed = pd.concat(df_list, keys=row_ix, axis=0)
         elif row_ix is None:
@@ -409,7 +422,7 @@ class VectorizedDF:
             X_mi_reconstructed = pd.concat(col_concats, keys=row_ix, axis=0)
 
         X_mi_index = X_mi_reconstructed.index
-        X_orig_row_index = self.X_multiindex.index
+        X_orig_row_index = self.X_mi_index
 
         flatten = col_multiindex == "flat" or (col_multiindex == "none" and force_flat)
         if flatten and isinstance(X_mi_reconstructed.columns, pd.MultiIndex):
diff --git a/sktime/forecasting/compose/_ensemble.py b/sktime/forecasting/compose/_ensemble.py
index f523479d6..e4211319b 100644
--- a/sktime/forecasting/compose/_ensemble.py
+++ b/sktime/forecasting/compose/_ensemble.py
@@ -364,8 +364,14 @@ class EnsembleForecaster(_HeterogenousEnsembleForecaster):
         """
         names, _ = self._check_forecasters()
         y_pred = pd.concat(self._predict_forecasters(fh, X), axis=1, keys=names)
-        y_pred = y_pred.groupby(level=1, axis=1).agg(
-            _aggregate, self.aggfunc, self.weights
+        y_pred = (
+            y_pred.T.groupby(level=1)
+            .agg(
+                lambda y, aggfunc, weights: _aggregate(y.T, aggfunc, weights),
+                self.aggfunc,
+                self.weights,
+            )
+            .T
         )
         return y_pred
 
diff --git a/sktime/forecasting/model_selection/_tune.py b/sktime/forecasting/model_selection/_tune.py
index 5c11825bf..05ff5c5be 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -401,10 +401,10 @@ class ForecastingGridSearchCV(BaseGridSearch):
 
     Parameters
     ----------
-    forecaster : estimator object
-        The estimator should implement the sktime or scikit-learn estimator
-        interface. Either the estimator must contain a "score" function,
-        or a scoring function must be passed.
+    forecaster : sktime forecaster, BaseForecaster instance or interface compatible
+        The forecaster to tune, must implement the sktime forecaster interface.
+        sklearn regressors can be used, but must first be converted to forecasters
+        via one of the reduction compositors, e.g., via ``make_reduction``
     cv : cross-validation generator or an iterable
         e.g. SlidingWindowSplitter()
     strategy : {"refit", "update", "no-update_params"}, optional, default="refit"
@@ -728,10 +728,10 @@ class ForecastingRandomizedSearchCV(BaseGridSearch):
 
     Parameters
     ----------
-    forecaster : estimator object
-        The estimator should implement the sktime or scikit-learn estimator
-        interface. Either the estimator must contain a "score" function,
-        or a scoring function must be passed.
+    forecaster : sktime forecaster, BaseForecaster instance or interface compatible
+        The forecaster to tune, must implement the sktime forecaster interface.
+        sklearn regressors can be used, but must first be converted to forecasters
+        via one of the reduction compositors, e.g., via ``make_reduction``
     cv : cross-validation generator or an iterable
         e.g. SlidingWindowSplitter()
     strategy : {"refit", "update", "no-update_params"}, optional, default="refit"
@@ -973,10 +973,10 @@ class ForecastingSkoptSearchCV(BaseGridSearch):
 
     Parameters
     ----------
-    forecaster : estimator object.
-        The estimator should implement the sktime or scikit-learn estimator interface.
-        Either the estimator must contain a "score" function,
-        or a scoring function must be passed.
+    forecaster : sktime forecaster, BaseForecaster instance or interface compatible
+        The forecaster to tune, must implement the sktime forecaster interface.
+        sklearn regressors can be used, but must first be converted to forecasters
+        via one of the reduction compositors, e.g., via ``make_reduction``
     cv : cross-validation generator or an iterable
         Splitter used for generating validation folds.
         e.g. SlidingWindowSplitter()
diff --git a/sktime/param_est/base.py b/sktime/param_est/base.py
index d1a430e76..9f4bb8335 100644
--- a/sktime/param_est/base.py
+++ b/sktime/param_est/base.py
@@ -381,23 +381,3 @@ class BaseParamFitter(BaseEstimator):
         # but looping to self.fit for now to avoid interface break
 
         return self
-
-    def _get_fitted_params(self):
-        """Get fitted parameters.
-
-        private _get_fitted_params, called from get_fitted_params
-
-        State required:
-            Requires state to be "fitted".
-
-        Returns
-        -------
-        fitted_params : dict
-        """
-        # default retrieves all self attributes ending in "_"
-        # and returns them with keys that have the "_" removed
-        fitted_params = [attr for attr in dir(self) if attr.endswith("_")]
-        fitted_params = [x for x in fitted_params if not x.startswith("_")]
-        fitted_param_dict = {p[:-1]: getattr(self, p) for p in fitted_params}
-
-        return fitted_param_dict
diff --git a/sktime/param_est/stationarity/_statsmodels.py b/sktime/param_est/stationarity/_statsmodels.py
index 21813a4b7..20bd57e64 100644
--- a/sktime/param_est/stationarity/_statsmodels.py
+++ b/sktime/param_est/stationarity/_statsmodels.py
@@ -158,7 +158,7 @@ class StationarityADF(BaseParamFitter):
 class StationarityKPSS(BaseParamFitter):
     """Test for stationarity via the Kwiatkowski-Phillips-Schmidt-Shin Test.
 
-    Uses ``statsmodels.tsa.stattools.kpss`` as a test for trend-stationairty,
+    Uses ``statsmodels.tsa.stattools.kpss`` as a test for trend-stationarity,
     and derives a boolean statement whether a series is (trend-)stationary.
 
     Also returns test results for the trend-stationarity test as fitted parameters.
diff --git a/sktime/performance_metrics/forecasting/probabilistic/_classes.py b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
index 2c6d07223..9208d73c3 100644
--- a/sktime/performance_metrics/forecasting/probabilistic/_classes.py
+++ b/sktime/performance_metrics/forecasting/probabilistic/_classes.py
@@ -119,16 +119,16 @@ class _BaseProbaForecastingErrorMetric(BaseForecastingErrorMetric):
 
         if isinstance(multioutput, str):
             if self.score_average and multioutput == "uniform_average":
-                out = float(out.mean(axis=1).iloc[0])  # average over all
+                out = out.mean(axis=1).iloc[0]  # average over all
             if self.score_average and multioutput == "raw_values":
-                out = out.groupby(axis=1, level=0).mean()  # average over scores
+                out = out.T.groupby(level=0).mean().T  # average over scores
             if not self.score_average and multioutput == "uniform_average":
-                out = out.groupby(axis=1, level=1).mean()  # average over variables
+                out = out.T.groupby(level=1).mean().T  # average over variables
             if not self.score_average and multioutput == "raw_values":
                 out = out  # don't average
         else:  # is np.array with weights
             if self.score_average:
-                out_raw = out.groupby(axis=1, level=0).mean()
+                out_raw = out.T.groupby(level=0).mean().T
                 out = out_raw.dot(multioutput)[0]
             else:
                 out = _groupby_dot(out, multioutput)
@@ -215,14 +215,14 @@ class _BaseProbaForecastingErrorMetric(BaseForecastingErrorMetric):
             if self.score_average and multioutput == "uniform_average":
                 out = out.mean(axis=1)  # average over all
             if self.score_average and multioutput == "raw_values":
-                out = out.groupby(axis=1, level=0).mean()  # average over scores
+                out = out.T.groupby(level=0).mean().T  # average over scores
             if not self.score_average and multioutput == "uniform_average":
-                out = out.groupby(axis=1, level=1).mean()  # average over variables
+                out = out.T.groupby(level=1).mean().T  # average over variables
             if not self.score_average and multioutput == "raw_values":
                 out = out  # don't average
         else:  # numpy array
             if self.score_average:
-                out_raw = out.groupby(axis=1, level=0).mean()
+                out_raw = out.T.groupby(level=0).mean().T
                 out = out_raw.dot(multioutput)
             else:
                 out = _groupby_dot(out, multioutput)
diff --git a/sktime/proba/base.py b/sktime/proba/base.py
index a1dd59f9a..8dc315094 100644
--- a/sktime/proba/base.py
+++ b/sktime/proba/base.py
@@ -29,6 +29,7 @@ class BaseDistribution(BaseObject):
         "approx_var_spl": 1000,  # sample size used in MC estimates of var
         "approx_energy_spl": 1000,  # sample size used in MC estimates of energy
         "approx_spl": 1000,  # sample size used in other MC estimates
+        "bisect_iter": 1000,  # max iters for bisection method in ppf
     }
 
     def __init__(self, index=None, columns=None):
@@ -42,15 +43,15 @@ class BaseDistribution(BaseObject):
     def loc(self):
         """Location indexer.
 
-        Use `my_distribution.loc[index]` for `pandas`-like row/column subsetting of
-        `BaseDistribution` descendants.
+        Use `my_distribution.loc[index]` for `pandas`-like row/column subsetting
+        of `BaseDistribution` descendants.
 
         `index` can be any `pandas` `loc` compatible index subsetter.
 
         `my_distribution.loc[index]` or `my_distribution.loc[row_index, col_index]`
-        subset `my_distribution` to rows defined by `row_index`, cols by `col_index`, to
-        exactly the same/cols rows as `pandas` `loc` would subset rows in
-        `my_distribution.index` and columns in `my_distribution.columns`.
+        subset `my_distribution` to rows defined by `row_index`, cols by `col_index`,
+        to exactly the same/cols rows as `pandas` `loc` would subset
+        rows in `my_distribution.index` and columns in `my_distribution.columns`.
         """
         return _Indexer(ref=self, method="_loc")
 
@@ -58,15 +59,15 @@ class BaseDistribution(BaseObject):
     def iloc(self):
         """Integer location indexer.
 
-        Use `my_distribution.iloc[index]` for `pandas`-like row/column subsetting of
-        `BaseDistribution` descendants.
+        Use `my_distribution.iloc[index]` for `pandas`-like row/column subsetting
+        of `BaseDistribution` descendants.
 
         `index` can be any `pandas` `iloc` compatible index subsetter.
 
         `my_distribution.iloc[index]` or `my_distribution.iloc[row_index, col_index]`
-        subset `my_distribution` to rows defined by `row_index`, cols by `col_index`, to
-        exactly the same/cols rows as `pandas` `iloc` would subset rows in
-        `my_distribution.index` and columns in `my_distribution.columns`.
+        subset `my_distribution` to rows defined by `row_index`, cols by `col_index`,
+        to exactly the same/cols rows as `pandas` `iloc` would subset
+        rows in `my_distribution.index` and columns in `my_distribution.columns`.
         """
         return _Indexer(ref=self, method="_iloc")
 
@@ -91,9 +92,12 @@ class BaseDistribution(BaseObject):
 
         subset_param_dict = {}
         for param, val in params.items():
-            arr = np.array(val)
-            if len(arr.shape) == 0:
-                subset_param_dict
+            if val is not None:
+                arr = np.array(val)
+            else:
+                arr = None
+            # if len(arr.shape) == 0:
+            # do nothing with arr
             if len(arr.shape) >= 1 and rowidx is not None:
                 arr = arr[rowidx]
             if len(arr.shape) >= 2 and colidx is not None:
@@ -283,10 +287,42 @@ class BaseDistribution(BaseObject):
         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."""
+        if self._has_implementation_of("cdf"):
+            from scipy.optimize import bisect
+
+            max_iter = self.get_tag("bisect_iter")
+            approx_method = (
+                "by using the bisection method (scipy.optimize.bisect) on "
+                f"the cdf, at {max_iter} maximum iterations"
+            )
+            warn(self._method_error_msg("cdf", fill_in=approx_method))
+
+            result = pd.DataFrame(index=p.index, columns=p.columns, dtype="float")
+            for ix in p.index:
+                for col in p.columns:
+                    d_ix = self.loc[[ix], [col]]
+                    p_ix = p.loc[ix, col]
+
+                    def opt_fun(x):
+                        """Optimization function, to find x s.t. cdf(x) = p_ix."""
+                        x = pd.DataFrame(x, index=[ix], columns=[col])  # noqa: B023
+                        return d_ix.cdf(x).values[0][0] - p_ix  # noqa: B023
+
+                    left_bd = -1e6
+                    right_bd = 1e6
+                    while opt_fun(left_bd) > 0:
+                        left_bd *= 10
+                    while opt_fun(right_bd) < 0:
+                        right_bd *= 10
+                    result.loc[ix, col] = bisect(
+                        opt_fun, left_bd, right_bd, maxiter=max_iter
+                    )
+            return result
+
         raise NotImplementedError(self._method_error_msg("ppf", "error"))
 
     def energy(self, x=None):
@@ -331,7 +367,8 @@ class BaseDistribution(BaseObject):
 
         # 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
 
@@ -354,7 +391,7 @@ class BaseDistribution(BaseObject):
         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.
@@ -377,10 +414,10 @@ class BaseDistribution(BaseObject):
         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.
+        r"""a-norm of pdf, defaults to 2-norm.
 
         computes a-norm of the entry marginal pdf, i.e.,
         :math:`\mathbb{E}[p_X(X)^{a-1}] = \int p(x)^a dx`,
@@ -409,7 +446,7 @@ class BaseDistribution(BaseObject):
 
         # 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)
@@ -431,7 +468,7 @@ class BaseDistribution(BaseObject):
         The `ppf` method also computes quantiles, but broadcasts differently, in
         `numpy` style closer to `tensorflow`.
         In contrast, this `quantile` method broadcasts
-        as forecaster `predict_quantiles`, i.e., columns first.
+        as ``sktime`` forecaster `predict_quantiles`, i.e., columns first.
 
         Parameters
         ----------
diff --git a/sktime/proba/empirical.py b/sktime/proba/empirical.py
index b19ab5f0e..f348c6e11 100644
--- a/sktime/proba/empirical.py
+++ b/sktime/proba/empirical.py
@@ -16,12 +16,12 @@ class Empirical(BaseDistribution):
     ----------
     spl : pd.DataFrame with pd.MultiIndex
         empirical sample
-        last (highest) index is time, first (lowest) index is sample
+        last (highest) index is instance, first (lowest) index is sample
     weights : pd.Series, with same index and length as spl, optional, default=None
         if not passed, ``spl`` is assumed to be unweighted
     time_indep : bool, optional, default=True
-        if True, ``sample`` will sample individual time indices independently
-        if False, ``sample`` will sample etire instances from ``spl``
+        if True, ``sample`` will sample individual instance indices independently
+        if False, ``sample`` will sample entire instances from ``spl``
     index : pd.Index, optional, default = RangeIndex
     columns : pd.Index, optional, default = RangeIndex
 
@@ -79,7 +79,6 @@ class Empirical(BaseDistribution):
 
         sorted = {}
         weights = {}
-        weights
         for t in times:
             sorted[t] = {}
             weights[t] = {}
@@ -124,8 +123,42 @@ class Empirical(BaseDistribution):
                     x_t = x.loc[ix, col]
                 else:
                     x_t = x
-                res.loc[ix, col] = func(spl=spl_t, weights=weights_t, x=x_t, **params)
-        return res.convert_dtypes()
+                res.at[ix, col] = func(spl=spl_t, weights=weights_t, x=x_t, **params)
+        return res.apply(pd.to_numeric)
+
+    def _iloc(self, rowidx=None, colidx=None):
+        index = self.index
+        columns = self.columns
+        weights = self.weights
+
+        spl_subset = self.spl
+
+        if rowidx is not None:
+            rowidx_loc = index[rowidx]
+            # subset multiindex to rowidx by last level
+            spl_subset = self.spl.loc[(slice(None), rowidx_loc), :]
+            if weights is not None:
+                weights_subset = weights.loc[(slice(None), rowidx_loc)]
+            else:
+                weights_subset = None
+            subs_rowidx = index[rowidx]
+        else:
+            subs_rowidx = index
+            weights_subset = weights
+
+        if colidx is not None:
+            spl_subset = spl_subset.iloc[:, colidx]
+            subs_colidx = columns[colidx]
+        else:
+            subs_colidx = columns
+
+        return Empirical(
+            spl_subset,
+            weights=weights_subset,
+            time_indep=self.time_indep,
+            index=subs_rowidx,
+            columns=subs_colidx,
+        )
 
     def energy(self, x=None):
         r"""Energy of self, w.r.t. self or a constant frame x.
@@ -162,9 +195,9 @@ class Empirical(BaseDistribution):
         """
         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)
@@ -186,11 +219,11 @@ class Empirical(BaseDistribution):
         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],
@@ -334,7 +367,7 @@ def _energy_np(spl, x=None, weights=None, assume_sorted=False):
 
     if x is None:
         cum_fwd = np.cumsum(weights[:-1])
-        cum_back = np.cumsum(weights[1::-1])[::-1]
+        cum_back = np.cumsum(weights[1:][::-1])[::-1]
         energy = 2 * np.sum(cum_fwd * cum_back * spl_diff)
     else:
         spl_diff = np.abs(spl - x)
diff --git a/sktime/proba/laplace.py b/sktime/proba/laplace.py
index 8f8345c61..274ce012b 100644
--- a/sktime/proba/laplace.py
+++ b/sktime/proba/laplace.py
@@ -43,7 +43,7 @@ class Laplace(BaseDistribution):
         # todo: untangle index handling
         # and broadcast of parameters.
         # move this functionality to the base class
-        self._mu, self._scale = self._get_bc_params()
+        self._mu, self._scale = self._get_bc_params(self.mu, self.scale)
         shape = self._mu.shape
 
         if index is None:
@@ -54,16 +54,6 @@ class Laplace(BaseDistribution):
 
         super().__init__(index=index, columns=columns)
 
-    def _get_bc_params(self):
-        """Fully broadcast parameters of self, given param shapes and index, columns."""
-        to_broadcast = [self.mu, self.scale]
-        if hasattr(self, "index") and self.index is not None:
-            to_broadcast += [self.index.to_numpy().reshape(-1, 1)]
-        if hasattr(self, "columns") and self.columns is not None:
-            to_broadcast += [self.columns.to_numpy()]
-        bc = np.broadcast_arrays(*to_broadcast)
-        return bc[0], bc[1]
-
     def energy(self, x=None):
         r"""Energy of self, w.r.t. self or a constant frame x.
 
diff --git a/sktime/proba/normal.py b/sktime/proba/normal.py
index 6de59c06e..3d51a19bb 100644
--- a/sktime/proba/normal.py
+++ b/sktime/proba/normal.py
@@ -44,7 +44,6 @@ class Normal(BaseDistribution):
         # todo: untangle index handling
         # and broadcast of parameters.
         # move this functionality to the base class
-        # 0.19.0?
         self._mu, self._sigma = self._get_bc_params(self.mu, self.sigma)
         shape = self._mu.shape
 
diff --git a/sktime/registry/_tags.py b/sktime/registry/_tags.py
index b40bb47b3..61da77999 100644
--- a/sktime/registry/_tags.py
+++ b/sktime/registry/_tags.py
@@ -492,6 +492,12 @@ ESTIMATOR_TAG_REGISTER = [
         "int",
         "sample size used in approximating other statistics if not available",
     ),
+    (
+        "bisect_iter",
+        "distribution",
+        "int",
+        "max iters for bisection method in ppf",
+    ),
     (
         "capability:multioutput",
         ["classifier", "regressor"],  # might need to add "early_classifier" here
diff --git a/sktime/transformations/panel/rocket/_multirocket.py b/sktime/transformations/panel/rocket/_multirocket.py
index 4cc186e17..3310eff46 100644
--- a/sktime/transformations/panel/rocket/_multirocket.py
+++ b/sktime/transformations/panel/rocket/_multirocket.py
@@ -8,6 +8,8 @@ import pandas as pd
 from sktime.datatypes import convert
 from sktime.transformations.base import BaseTransformer
 
+__author__ = ["ChangWeiTan", "fstinner", "angus924"]
+
 
 class MultiRocket(BaseTransformer):
     """Multi RandOm Convolutional KErnel Transform (MultiRocket).
diff --git a/sktime/transformations/panel/rocket/_multirocket_multivariate.py b/sktime/transformations/panel/rocket/_multirocket_multivariate.py
index 1fe4813f4..dbc86b338 100644
--- a/sktime/transformations/panel/rocket/_multirocket_multivariate.py
+++ b/sktime/transformations/panel/rocket/_multirocket_multivariate.py
@@ -5,6 +5,8 @@ import pandas as pd
 
 from sktime.transformations.base import BaseTransformer
 
+__author__ = ["ChangWeiTan", "fstinner", "angus924"]
+
 
 class MultiRocketMultivariate(BaseTransformer):
     """Multi RandOm Convolutional KErnel Transform (MultiRocket).
diff --git a/sktime/utils/index.py b/sktime/utils/index.py
new file mode 100644
index 000000000..c9c4cde9f
--- /dev/null
+++ b/sktime/utils/index.py
@@ -0,0 +1,20 @@
+"""Utility functions for working with indices."""
+
+import numpy as np
+
+
+def random_ss_ix(ix, size, replace=True):
+    """Randomly uniformly sample indices from a list of indices.
+
+    Parameters
+    ----------
+    ix : pd.Index or subsettable iterable via getitem
+        list of indices to sample from
+    size : int
+        number of indices to sample
+    replace : bool, default=True
+        whether to sample with replacement
+    """
+    a = range(len(ix))
+    ixs = ix[np.random.choice(a, size=size, replace=replace)]
+    return ixs
diff --git a/sktime/utils/parallel.py b/sktime/utils/parallel.py
index 6a8fb4090..fbfb9de63 100644
--- a/sktime/utils/parallel.py
+++ b/sktime/utils/parallel.py
@@ -108,13 +108,13 @@ def _parallelize_joblib(fun, iter, meta, backend, backend_params):
             raise ValueError(
                 '"joblib" was selected as first layer parallelization backend, '
                 "but no backend string was "
-                '"passed in the backend parameters dict, e.g., "spark". '
+                'passed in the backend parameters dict, e.g., "spark". '
                 "Please specify a backend to joblib as a key-value pair "
                 "in the backend_params arg or the backend:parallel:params config "
                 'when using "joblib". '
                 'For clarity, "joblib" should only be used for two-layer '
                 "backend dispatch, where the first layer is joblib, "
-                "and the second layer is a custom backend of joblib, e.g., spark."
+                "and the second layer is a custom backend of joblib, e.g., spark. "
                 "For first-party joblib backends, please use the backend string "
                 'of sktime directly, e.g., by specifying "multiprocessing" or "loky".'
             )
@@ -149,11 +149,29 @@ def _parallelize_dask(fun, iter, meta, backend, backend_params):
 para_dict["dask"] = _parallelize_dask
 
 
-def _get_parallel_test_fixtures():
+def _get_parallel_test_fixtures(naming="estimator"):
     """Return fixtures for parallelization tests.
 
     Returns a list of parameter fixtures, where each fixture
     is a dict with keys "backend" and "backend_params".
+
+    Parameters
+    ----------
+    naming : str, optional
+        naming convention for the parameters, one of
+
+        "estimator": for use in estimator constructors,
+        ``backend`` and ``backend_params``
+        "config": for use in ``set_config``,
+        ``backend:parallel`` and ``backend:parallel:params``
+
+    Returns
+    -------
+    fixtures : list of dict
+        list of backend parameter fixtures
+        keys depend on ``naming`` parameter, see above
+        values are backend strings and backend parameter dicts
+        only backends that are available in the environment are included
     """
     from sktime.utils.validation._dependencies import _check_soft_dependencies
 


Test patch:
diff --git a/sktime/classification/tests/test_all_classifiers.py b/sktime/classification/tests/test_all_classifiers.py
index cda9a0245..173edc08a 100644
--- a/sktime/classification/tests/test_all_classifiers.py
+++ b/sktime/classification/tests/test_all_classifiers.py
@@ -229,7 +229,8 @@ class TestAllClassifiers(ClassifierFixtureGenerator, QuickTester):
         assert isinstance(y_pred, pd.DataFrame)
         assert y_pred.shape == y_mult.shape
 
-        vectorized = estimator_instance.get_tag("capability:multioutput")
+        # the estimator vectorizes iff it does not have the multioutput capability
+        vectorized = not estimator_instance.get_tag("capability:multioutput")
         if vectorized:
             assert hasattr(estimator_instance, "classifiers_")
             assert isinstance(estimator_instance.classifiers_, pd.DataFrame)
diff --git a/sktime/forecasting/base/tests/test_base.py b/sktime/forecasting/base/tests/test_base.py
index b2dc5c0d1..a0d04b69b 100644
--- a/sktime/forecasting/base/tests/test_base.py
+++ b/sktime/forecasting/base/tests/test_base.py
@@ -20,6 +20,7 @@ from sktime.forecasting.var import VAR
 from sktime.utils._testing.hierarchical import _make_hierarchical
 from sktime.utils._testing.panel import _make_panel
 from sktime.utils._testing.series import _make_series
+from sktime.utils.parallel import _get_parallel_test_fixtures
 from sktime.utils.validation._dependencies import _check_estimator_deps
 
 PANEL_MTYPES = ["pd-multiindex", "nested_univ", "numpy3D"]
@@ -30,7 +31,7 @@ HIER_MTYPES = ["pd_multiindex_hier"]
     not _check_estimator_deps(ARIMA, severity="none"),
     reason="skip test if required soft dependency for ARIMA not available",
 )
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", _get_parallel_test_fixtures("config"))
 @pytest.mark.parametrize("mtype", PANEL_MTYPES)
 def test_vectorization_series_to_panel(mtype, backend):
     """Test that forecaster vectorization works for Panel data.
@@ -43,7 +44,7 @@ def test_vectorization_series_to_panel(mtype, backend):
     y = _make_panel(n_instances=n_instances, random_state=42, return_mtype=mtype)
 
     f = ARIMA()
-    f.set_config(**{"backend:parallel": backend})
+    f.set_config(**backend.copy())
     y_pred = f.fit(y).predict([1, 2, 3])
     valid, _, metadata = check_is_mtype(
         y_pred, mtype, return_metadata=True, msg_return_dict="list"
diff --git a/sktime/forecasting/model_selection/tests/test_tune.py b/sktime/forecasting/model_selection/tests/test_tune.py
index fa52904df..e1cb65dd6 100644
--- a/sktime/forecasting/model_selection/tests/test_tune.py
+++ b/sktime/forecasting/model_selection/tests/test_tune.py
@@ -325,7 +325,7 @@ def test_skoptcv_multiple_forecaster():
     assert len(sscv.cv_results_) == 5
 
 
-BACKEND_TEST = _get_parallel_test_fixtures()
+BACKEND_TEST = _get_parallel_test_fixtures("estimator")
 
 
 @pytest.mark.skipif(
diff --git a/sktime/proba/tests/test_all_distrs.py b/sktime/proba/tests/test_all_distrs.py
index ae5a3334c..0a12fc98c 100644
--- a/sktime/proba/tests/test_all_distrs.py
+++ b/sktime/proba/tests/test_all_distrs.py
@@ -1,7 +1,7 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Tests for BaseDistribution API points."""
 
-__author__ = ["fkiraly"]
+__author__ = ["fkiraly", "Alex-JG3"]
 
 import numpy as np
 import pandas as pd
@@ -9,6 +9,7 @@ import pytest
 
 from sktime.datatypes import check_is_mtype
 from sktime.tests.test_all_estimators import BaseFixtureGenerator, QuickTester
+from sktime.utils.index import random_ss_ix
 
 
 class DistributionFixtureGenerator(BaseFixtureGenerator):
@@ -54,12 +55,16 @@ METHODS_ROWWISE = ["energy"]  # results in one column
 
 
 class TestAllDistributions(DistributionFixtureGenerator, QuickTester):
-    """Module level tests for all sktime parameter fitters."""
+    """Module level tests for all sktime distributions."""
 
-    def test_sample(self, estimator_instance):
+    @pytest.mark.parametrize("shuffled", [False, True])
+    def test_sample(self, estimator_instance, shuffled):
         """Test sample expected return."""
         d = estimator_instance
 
+        if shuffled:
+            d = _shuffle_distr(d)
+
         res = d.sample()
 
         assert d.shape == res.shape
@@ -72,36 +77,50 @@ class TestAllDistributions(DistributionFixtureGenerator, QuickTester):
         assert (res_panel.index == dummy_panel.index).all()
         assert (res_panel.columns == dummy_panel.columns).all()
 
+    @pytest.mark.parametrize("shuffled", [False, True])
     @pytest.mark.parametrize("method", METHODS_SCALAR, ids=METHODS_SCALAR)
-    def test_methods_scalar(self, estimator_instance, method):
+    def test_methods_scalar(self, estimator_instance, method, shuffled):
         """Test expected return of scalar methods."""
         if not _has_capability(estimator_instance, method):
             return None
 
         d = estimator_instance
-        res = getattr(estimator_instance, method)()
+        if shuffled:
+            d = _shuffle_distr(d)
+
+        res = getattr(d, method)()
 
         _check_output_format(res, d, method)
 
+    @pytest.mark.parametrize("shuffled", [False, True])
     @pytest.mark.parametrize("method", METHODS_X, ids=METHODS_X)
-    def test_methods_x(self, estimator_instance, method):
+    def test_methods_x(self, estimator_instance, method, shuffled):
         """Test expected return of methods that take sample-like argument."""
         if not _has_capability(estimator_instance, method):
             return None
 
         d = estimator_instance
+
+        if shuffled:
+            d = _shuffle_distr(d)
+
         x = d.sample()
-        res = getattr(estimator_instance, method)(x)
+        res = getattr(d, method)(x)
 
         _check_output_format(res, d, method)
 
+    @pytest.mark.parametrize("shuffled", [False, True])
     @pytest.mark.parametrize("method", METHODS_P, ids=METHODS_P)
-    def test_methods_p(self, estimator_instance, method):
+    def test_methods_p(self, estimator_instance, method, shuffled):
         """Test expected return of methods that take percentage-like argument."""
         if not _has_capability(estimator_instance, method):
             return None
 
         d = estimator_instance
+
+        if shuffled:
+            d = _shuffle_distr(d)
+
         np_unif = np.random.uniform(size=d.shape)
         p = pd.DataFrame(np_unif, index=d.index, columns=d.columns)
         res = getattr(estimator_instance, method)(p)
@@ -130,6 +149,63 @@ class TestAllDistributions(DistributionFixtureGenerator, QuickTester):
         res = d.quantile(q)
         _check_quantile_output(res, q)
 
+    @pytest.mark.parametrize("subset_row", [True, False])
+    @pytest.mark.parametrize("subset_col", [True, False])
+    def test_subsetting(self, estimator_instance, subset_row, subset_col):
+        """Test subsetting of distribution."""
+        d = estimator_instance
+
+        if subset_row:
+            ix_loc = random_ss_ix(d.index, 3)
+            ix_iloc = d.index.get_indexer(ix_loc)
+        else:
+            ix_loc = d.index
+            ix_iloc = pd.RangeIndex(len(d.index))
+
+        if subset_col:
+            iy_loc = random_ss_ix(d.columns, 1)
+            iy_iloc = d.columns.get_indexer(iy_loc)
+        else:
+            iy_loc = d.columns
+            iy_iloc = pd.RangeIndex(len(d.columns))
+
+        res_loc = d.loc[ix_loc, iy_loc]
+
+        assert isinstance(res_loc, type(d))
+        assert res_loc.shape == (len(ix_loc), len(iy_loc))
+        assert (res_loc.index == ix_loc).all()
+        assert (res_loc.columns == iy_loc).all()
+
+        res_iloc = d.iloc[ix_iloc, iy_iloc]
+
+        assert isinstance(res_iloc, type(d))
+        assert res_iloc.shape == (len(ix_iloc), len(iy_iloc))
+        assert (res_iloc.index == ix_loc).all()
+        assert (res_iloc.columns == iy_loc).all()
+
+    def test_log_pdf_and_pdf(self, estimator_instance):
+        """Test that the log of the pdf and log_pdf function are similar."""
+        d = estimator_instance
+        capabilities_exact = d.get_tags()["capabilities:exact"]
+
+        if "log_pdf" not in capabilities_exact or "pdf" not in capabilities_exact:
+            return
+        x = d.sample()
+        pdf = d.pdf(x)
+        log_pdf = d.log_pdf(x)
+        assert np.allclose(np.log(pdf), log_pdf)
+
+    def test_ppf_and_cdf(self, estimator_instance):
+        """Test that the ppf is the inverse of the cdf."""
+        d = estimator_instance
+        capabilities_exact = d.get_tags()["capabilities:exact"]
+
+        if "ppf" not in capabilities_exact or "cdf" not in capabilities_exact:
+            return
+        x = d.sample()
+        x_approx = d.ppf(d.cdf(x))
+        assert np.allclose(x.values, x_approx.values)
+
 
 def _check_output_format(res, dist, method):
     """Check output format expectations for BaseDistribution tests."""
@@ -144,3 +220,16 @@ def _check_output_format(res, dist, method):
 
     if method in METHODS_SCALAR_POS or method in METHODS_X_POS:
         assert (res >= 0).all().all()
+
+    if isinstance(res, pd.DataFrame):
+        assert res.apply(pd.api.types.is_numeric_dtype).all()
+    elif isinstance(res, pd.Series):
+        assert pd.api.types.is_numeric_dtype(res)
+    else:
+        raise TypeError("res must be a pandas DataFrame or Series.")
+
+
+def _shuffle_distr(d):
+    """Shuffle distribution row index."""
+    shuffled_index = pd.DataFrame(d.index).sample(frac=1).index
+    return d.loc[shuffled_index]
diff --git a/sktime/transformations/tests/test_base.py b/sktime/transformations/tests/test_base.py
index b18b1e918..be3b7f85b 100644
--- a/sktime/transformations/tests/test_base.py
+++ b/sktime/transformations/tests/test_base.py
@@ -35,15 +35,14 @@ from sktime.utils._testing.scenarios_transformers import (
     TransformerFitTransformSeriesUnivariate,
 )
 from sktime.utils._testing.series import _make_series
+from sktime.utils.parallel import _get_parallel_test_fixtures
 from sktime.utils.validation._dependencies import _check_soft_dependencies
 
 # other scenarios that might be needed later in development:
 # TransformerFitTransformPanelUnivariateWithClassY,
 
 # list of parallelization backends to test
-BACKENDS = [None, "multiprocessing", "loky", "threading"]
-if _check_soft_dependencies("dask", severity="none"):
-    BACKENDS.append("dask")
+BACKENDS = _get_parallel_test_fixtures("config")
 
 
 def inner_X_scitypes(est):
@@ -204,7 +203,7 @@ def test_panel_in_panel_out_not_supported_but_series(backend):
     # one example for a transformer which supports Series internally but not Panel
     cls = BoxCoxTransformer
     est = cls.create_test_instance()
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing more scitypes)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -274,7 +273,7 @@ def test_panel_in_primitives_out_not_supported_fit_in_transform(backend):
     # one example for a transformer which supports Series internally but not Panel
     cls = SummaryTransformer
     est = cls.create_test_instance()
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing more scitypes)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -391,7 +390,7 @@ def test_hierarchical_in_hierarchical_out_not_supported_but_series(backend):
     # one example for a transformer which supports Series internally
     cls = BoxCoxTransformer
     est = cls.create_test_instance()
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing more scitypes)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -468,7 +467,7 @@ def test_vectorization_multivariate_no_row_vectorization(backend):
     # one example for a transformer which supports Series internally
     cls = BoxCoxTransformer
     est = cls.create_test_instance()
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing multivariate functionality)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -509,7 +508,7 @@ def test_vectorization_multivariate_and_hierarchical(backend):
     # one example for a transformer which supports Series internally
     cls = BoxCoxTransformer
     est = cls.create_test_instance()
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing more scitypes)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -552,7 +551,7 @@ def test_vectorization_multivariate_no_row_vectorization_empty_fit(backend):
     # one example for a transformer which supports Series internally
     cls = BoxCoxTransformer
     est = FitInTransform(cls.create_test_instance())
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing multivariate functionality)
     #   (then this is not a failure of cls, but we need to choose another example)
@@ -593,7 +592,7 @@ def test_vectorization_multivariate_and_hierarchical_empty_fit(backend):
     # one example for a transformer which supports Series internally
     cls = BoxCoxTransformer
     est = FitInTransform(cls.create_test_instance())
-    est.set_config(**{"backend:parallel": backend})
+    est.set_config(**backend.copy())
     # ensure cls is a good example, if this fails, choose another example
     #   (if this changes, it may be due to implementing more scitypes)
     #   (then this is not a failure of cls, but we need to choose another example)
diff --git a/sktime/utils/_testing/deep_equals.py b/sktime/utils/_testing/deep_equals.py
index 0616cc81a..c114c5d19 100644
--- a/sktime/utils/_testing/deep_equals.py
+++ b/sktime/utils/_testing/deep_equals.py
@@ -180,7 +180,10 @@ def deep_equals(x, y, return_msg=False):
     elif isinstance(x, np.ndarray):
         if x.dtype != y.dtype:
             return ret(False, f".dtype, x.dtype = {x.dtype} != y.dtype = {y.dtype}")
-        return ret(np.array_equal(x, y, equal_nan=True), ".values")
+        if x.dtype in ["object", "str"]:
+            return ret(np.array_equal(x, y), ".values")
+        else:
+            return ret(np.array_equal(x, y, equal_nan=True), ".values")
     # recursion through lists, tuples and dicts
     elif isinstance(x, (list, tuple)):
         return ret(*_tuple_equals(x, y, return_msg=True))
diff --git a/sktime/utils/_testing/tests/test_deep_equals.py b/sktime/utils/_testing/tests/test_deep_equals.py
index 9d16330eb..16be1e9c2 100644
--- a/sktime/utils/_testing/tests/test_deep_equals.py
+++ b/sktime/utils/_testing/tests/test_deep_equals.py
@@ -26,6 +26,10 @@ EXAMPLES = [
     ForecastingHorizon([1, 2, 3], is_relative=False),
     {"foo": [42], "bar": pd.Series([1, 2])},
     {"bar": [42], "foo": pd.Series([1, 2])},
+    pd.Index([1, 2, 3]),
+    pd.Index([2, 3, 4]),
+    np.array([0.1, 1], dtype="object"),
+    np.array([0.2, 1], dtype="object"),
 ]