Repository: sktime/sktime
Pull request #5711 (2024-01-07 16:58:57)

[DOC] improve clarity in describing `strategy="refit"` in forecasting tuners' docstrings.
This PR improves clarity in describing `strategy="refit"` for forecasting tuners.

Inspired by discussion here: https://github.com/sktime/sktime/discussions/5704


Patch:
diff --git a/README.md b/README.md
index a435ea1b1..bd4cae0c5 100644
--- a/README.md
+++ b/README.md
@@ -66,8 +66,6 @@ Our objective is to enhance the interoperability and usability of the time serie
 
 sktime also provides **interfaces to related libraries**, for example [scikit-learn], [statsmodels], [tsfresh], [PyOD], and [fbprophet], among others.
 
-For **deep learning**, see our companion package: [sktime-dl](https://github.com/sktime/sktime-dl).
-
 [statsmodels]: https://www.statsmodels.org/stable/index.html
 [tsfresh]: https://tsfresh.readthedocs.io/en/latest/
 [pyod]: https://pyod.readthedocs.io/en/latest/
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/forecasting.rst b/docs/source/api_reference/forecasting.rst
index 116f4eade..e7b87652b 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -72,6 +72,7 @@ Use ``make_reduction`` for easy specification.
 
     DirectTabularRegressionForecaster
     DirectTimeSeriesRegressionForecaster
+    DirectReductionForecaster
     MultioutputTabularRegressionForecaster
     MultioutputTimeSeriesRegressionForecaster
     RecursiveTabularRegressionForecaster
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/docs/source/api_reference/transformations.rst b/docs/source/api_reference/transformations.rst
index 0bbc48474..28028133a 100644
--- a/docs/source/api_reference/transformations.rst
+++ b/docs/source/api_reference/transformations.rst
@@ -242,6 +242,7 @@ Lagging
     :template: class.rst
 
     Lag
+    ReducerTransform
 
 Element-wise transforms
 ~~~~~~~~~~~~~~~~~~~~~~~
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..0a68e35b9 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -401,16 +401,16 @@ 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"
         data ingestion strategy in fitting cv, passed to `evaluate` internally
         defines the ingestion mode when the forecaster sees new data when window expands
-        "refit" = forecaster is refitted to each training window
+        "refit" = a new copy of the forecaster is fitted to each training window
         "update" = forecaster is updated with training window data, in sequence provided
         "no-update_params" = fit to first training window, re-used without fit or update
     update_behaviour : str, optional, default = "full_refit"
@@ -728,16 +728,16 @@ 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"
         data ingestion strategy in fitting cv, passed to `evaluate` internally
         defines the ingestion mode when the forecaster sees new data when window expands
-        "refit" = forecaster is refitted to each training window
+        "refit" = a new copy of the forecaster is fitted to each training window
         "update" = forecaster is updated with training window data, in sequence provided
         "no-update_params" = fit to first training window, re-used without fit or update
     update_behaviour: str, optional, default = "full_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()
@@ -1034,7 +1034,7 @@ class ForecastingSkoptSearchCV(BaseGridSearch):
     strategy : {"refit", "update", "no-update_params"}, optional, default="refit"
         data ingestion strategy in fitting cv, passed to `evaluate` internally
         defines the ingestion mode when the forecaster sees new data when window expands
-        "refit" = forecaster is refitted to each training window
+        "refit" = a new copy of the forecaster is fitted to each training window
         "update" = forecaster is updated with training window data, in sequence provided
         "no-update_params" = fit to first training window, re-used without fit or update
     update_behaviour: str, optional, default = "full_refit"
diff --git a/sktime/forecasting/tbats.py b/sktime/forecasting/tbats.py
index 9477bf76b..a346bacd3 100644
--- a/sktime/forecasting/tbats.py
+++ b/sktime/forecasting/tbats.py
@@ -134,7 +134,7 @@ class TBATS(_TbatsAdapter):
         -------
         params : dict or list of dict
         """
-        params = {
+        params1 = {
             "use_box_cox": False,
             "use_trend": False,
             "use_damped_trend": False,
@@ -142,4 +142,12 @@ class TBATS(_TbatsAdapter):
             "use_arma_errors": False,
             "n_jobs": 1,
         }
-        return params
+        params2 = {
+            "use_box_cox": False,
+            "use_trend": True,
+            "use_damped_trend": True,
+            "sp": [],
+            "use_arma_errors": True,
+            "n_jobs": 2,
+        }
+        return [params1, params2]
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 4e09393af..be3b7f85b 100644
--- a/sktime/transformations/tests/test_base.py
+++ b/sktime/transformations/tests/test_base.py
@@ -35,11 +35,15 @@ 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 = _get_parallel_test_fixtures("config")
+
 
 def inner_X_scitypes(est):
     """Return list of scitypes supported by class est, as list of str."""
@@ -183,7 +187,7 @@ def test_panel_in_panel_out_supported():
     # todo: possibly, add mtype check, use metadata return
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_panel_in_panel_out_not_supported_but_series(backend):
     """Test that fit/transform runs and returns the correct output type.
 
@@ -199,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)
@@ -253,7 +257,7 @@ def test_series_in_primitives_out_supported_fit_in_transform():
     assert len(Xt) == 1
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_panel_in_primitives_out_not_supported_fit_in_transform(backend):
     """Test that fit/transform runs and returns the correct output type.
 
@@ -269,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)
@@ -370,7 +374,7 @@ def test_panel_in_primitives_out_supported_with_y_in_fit_but_not_transform():
     assert len(Xt) == 7
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_hierarchical_in_hierarchical_out_not_supported_but_series(backend):
     """Test that fit/transform runs and returns the correct output type.
 
@@ -386,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)
@@ -444,7 +448,7 @@ def test_hierarchical_in_hierarchical_out_not_supported_but_series_fit_in_transf
     assert len(Xt) == 2 * 4 * 12
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_vectorization_multivariate_no_row_vectorization(backend):
     """Test that multivariate vectorization of univariate transformers works.
 
@@ -463,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)
@@ -485,7 +489,7 @@ def test_vectorization_multivariate_no_row_vectorization(backend):
     assert len(Xt.columns) == len(scenario.args["fit"]["X"].columns)
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_vectorization_multivariate_and_hierarchical(backend):
     """Test that fit/transform runs and returns the correct output type.
 
@@ -504,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)
@@ -528,7 +532,7 @@ def test_vectorization_multivariate_and_hierarchical(backend):
     assert len(Xt.columns) == len(scenario.args["fit"]["X"].columns)
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_vectorization_multivariate_no_row_vectorization_empty_fit(backend):
     """Test that multivariate vectorization of univariate transformers works.
 
@@ -547,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)
@@ -569,7 +573,7 @@ def test_vectorization_multivariate_no_row_vectorization_empty_fit(backend):
     assert len(Xt.columns) == len(scenario.args["fit"]["X"].columns)
 
 
-@pytest.mark.parametrize("backend", [None, "joblib", "loky", "threading"])
+@pytest.mark.parametrize("backend", BACKENDS)
 def test_vectorization_multivariate_and_hierarchical_empty_fit(backend):
     """Test that fit/transform runs and returns the correct output type.
 
@@ -588,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"),
 ]