Repository: sktime/sktime
Pull request #5726 (2024-01-11 00:43:26)

[BUG] fix `ForecastingHorizon.get_expected_pred_idx` `sort_time`
This fixes an unreported bug where the `sort_time` argument of `ForecastingHorizon.get_expected_pred_idx` would behave erratically if th original time series was not sorted.


Patch:
diff --git a/docs/source/api_reference/classification.rst b/docs/source/api_reference/classification.rst
index 9a2d3ee19..b003e5097 100644
--- a/docs/source/api_reference/classification.rst
+++ b/docs/source/api_reference/classification.rst
@@ -23,6 +23,17 @@ Composition
     SklearnClassifierPipeline
     MultiplexClassifier
 
+Model selection and tuning
+--------------------------
+
+.. currentmodule:: sktime.classification.model_selection
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: class.rst
+
+    TSCGridSearchCV
+
 Ensembles
 ---------
 
@@ -182,8 +193,15 @@ Shapelet-based
     MrSEQL
     MrSQM
 
-sklearn
--------
+
+sklearn classifiers
+-------------------
+
+This section contains classifiers which are not time series classifiers but
+simple tabular classifiers in ``sklearn`` compatible API.
+
+They are used internally in time series classifiers, but can also be used
+directly in a tabular setting.
 
 .. currentmodule:: sktime.classification.sklearn
 
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/classification/model_selection/__init__.py b/sktime/classification/model_selection/__init__.py
new file mode 100644
index 000000000..72343b457
--- /dev/null
+++ b/sktime/classification/model_selection/__init__.py
@@ -0,0 +1,6 @@
+"""Tuning of time series classifiers."""
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+
+from sktime.classification.model_selection._tune import TSCGridSearchCV
+
+__all__ = ["TSCGridSearchCV"]
diff --git a/sktime/classification/model_selection/_tune.py b/sktime/classification/model_selection/_tune.py
new file mode 100644
index 000000000..add7ea31b
--- /dev/null
+++ b/sktime/classification/model_selection/_tune.py
@@ -0,0 +1,442 @@
+"""Tuning for time series classifiers."""
+
+__author__ = ["fkiraly", "achieveordie"]
+
+import numpy as np
+from sklearn.model_selection import GridSearchCV
+
+from sktime.classification._delegate import _DelegatedClassifier
+
+
+class TSCGridSearchCV(_DelegatedClassifier):
+    """Exhaustive search over specified parameter values for an estimator.
+
+    Adapts sklearn GridSearchCV for sktime time series classifiers
+
+    Optimizes hyper-parameters of `estimators` by exhaustive grid search.
+
+    Parameters
+    ----------
+    estimator : estimator object
+        This is assumed to implement the scikit-learn estimator interface.
+        Either estimator needs to provide a ``score`` function,
+        or ``scoring`` must be passed.
+
+    param_grid : dict or list of dictionaries
+        Dictionary with parameters names (`str`) as keys and lists of
+        parameter settings to try as values, or a list of such
+        dictionaries, in which case the grids spanned by each dictionary
+        in the list are explored. This enables searching over any sequence
+        of parameter settings.
+
+    scoring : str, callable, list, tuple or dict, default=None
+        Strategy to evaluate the performance of the cross-validated model on
+        the test set.
+
+        If `scoring` represents a single score, one can use:
+
+        - a single string (see :ref:`scoring_parameter`);
+        - a callable (see :ref:`scoring`) that returns a single value.
+
+        If `scoring` represents multiple scores, one can use:
+
+        - a list or tuple of unique strings;
+        - a callable returning a dictionary where the keys are the metric
+          names and the values are the metric scores;
+        - a dictionary with metric names as keys and callables a values.
+
+    n_jobs : int, default=None
+        Number of jobs to run in parallel.
+        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
+        ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
+        for more details.
+
+    refit : bool, str, or callable, default=True
+        Refit an estimator using the best found parameters on the whole
+        dataset. If ``False``, the ``predict`` and ``predict_probab`` will not work.
+
+        For multiple metric evaluation, this needs to be a `str` denoting the
+        scorer that would be used to find the best parameters for refitting
+        the estimator at the end.
+
+        Where there are considerations other than maximum score in
+        choosing a best estimator, ``refit`` can be set to a function which
+        returns the selected ``best_index_`` given ``cv_results_``. In that
+        case, the ``best_estimator_`` and ``best_params_`` will be set
+        according to the returned ``best_index_`` while the ``best_score_``
+        attribute will not be available.
+
+        The refitted estimator is made available at the ``best_estimator_``
+        attribute and permits using ``predict`` directly on this
+        ``GridSearchCV`` instance.
+
+        Also for multiple metric evaluation, the attributes ``best_index_``,
+        ``best_score_`` and ``best_params_`` will only be available if
+        ``refit`` is set and all of them will be determined w.r.t this specific
+        scorer.
+
+        See ``scoring`` parameter to know more about multiple metric
+        evaluation.
+
+    cv : int, cross-validation generator or an iterable, default=None
+        Determines the cross-validation splitting strategy.
+        Possible inputs for cv are:
+
+        - None, to use the default 5-fold cross validation,
+        - integer, to specify the number of folds in a `(Stratified)KFold`,
+        - :term:`CV splitter`,
+        - An iterable yielding (train, test) splits as arrays of indices.
+
+        For integer/None inputs, if the estimator is a classifier and ``y`` is
+        either binary or multiclass, :class:`StratifiedKFold` is used. In all
+        other cases, :class:`KFold` is used. These splitters are instantiated
+        with `shuffle=False` so the splits will be the same across calls.
+
+        Refer :ref:`User Guide <cross_validation>` for the various
+        cross-validation strategies that can be used here.
+
+    verbose : int
+        Controls the verbosity: the higher, the more messages.
+
+        - >1 : the computation time for each fold and parameter candidate is
+          displayed;
+        - >2 : the score is also displayed;
+        - >3 : the fold and candidate parameter indexes are also displayed
+          together with the starting time of the computation.
+
+    pre_dispatch : int, or str, default='2*n_jobs'
+        Controls the number of jobs that get dispatched during parallel
+        execution. Reducing this number can be useful to avoid an
+        explosion of memory consumption when more jobs get dispatched
+        than CPUs can process. This parameter can be:
+
+            - None, in which case all the jobs are immediately
+              created and spawned. Use this for lightweight and
+              fast-running jobs, to avoid delays due to on-demand
+              spawning of the jobs
+
+            - An int, giving the exact number of total jobs that are
+              spawned
+
+            - A str, giving an expression as a function of n_jobs,
+              as in '2*n_jobs'
+
+    error_score : 'raise' or numeric, default=np.nan
+        Value to assign to the score if an error occurs in estimator fitting.
+        If set to 'raise', the error is raised. If a numeric value is given,
+        FitFailedWarning is raised. This parameter does not affect the refit
+        step, which will always raise the error.
+
+    return_train_score : bool, default=False
+        If ``False``, the ``cv_results_`` attribute will not include training
+        scores.
+        Computing training scores is used to get insights on how different
+        parameter settings impact the overfitting/underfitting trade-off.
+        However computing the scores on the training set can be computationally
+        expensive and is not strictly required to select the parameters that
+        yield the best generalization performance.
+
+    tune_by_variable : bool, optional (default=False)
+        Whether to tune parameter by each time series variable separately,
+        in case of multivariate data passed to the tuning estimator.
+        Only applies if time series passed are strictly multivariate.
+        If True, clones of the estimator will be fit to each variable separately,
+        and are available in fields of the classifiers_ attribute.
+        Has the same effect as applying ColumnEnsembleClassifier wrapper to self.
+        If False, the same best parameter is selected for all variables.
+
+    Attributes
+    ----------
+    cv_results_ : dict of numpy (masked) ndarrays
+        A dict with keys as column headers and values as columns, that can be
+        imported into a pandas ``DataFrame``.
+
+        For multi-metric evaluation, the scores for all the scorers are
+        available in the ``cv_results_`` dict at the keys ending with that
+        scorer's name (``'_<scorer_name>'``) instead of ``'_score'`` shown
+        above. ('split0_test_precision', 'mean_train_precision' etc.)
+
+    best_estimator_ : estimator
+        Estimator that was chosen by the search, i.e. estimator
+        which gave highest score (or smallest loss if specified)
+        on the left out data. Not available if ``refit=False``.
+
+        See ``refit`` parameter for more information on allowed values.
+
+    best_score_ : float
+        Mean cross-validated score of the best_estimator
+
+        For multi-metric evaluation, this is present only if ``refit`` is
+        specified.
+
+        This attribute is not available if ``refit`` is a function.
+
+    best_params_ : dict
+        Parameter setting that gave the best results on the hold out data.
+
+        For multi-metric evaluation, this is present only if ``refit`` is
+        specified.
+
+    best_index_ : int
+        The index (of the ``cv_results_`` arrays) which corresponds to the best
+        candidate parameter setting.
+
+        The dict at ``search.cv_results_['params'][search.best_index_]`` gives
+        the parameter setting for the best model, that gives the highest
+        mean score (``search.best_score_``).
+
+        For multi-metric evaluation, this is present only if ``refit`` is
+        specified.
+
+    scorer_ : function or a dict
+        Scorer function used on the held out data to choose the best
+        parameters for the model.
+
+        For multi-metric evaluation, this attribute holds the validated
+        ``scoring`` dict which maps the scorer key to the scorer callable.
+
+    n_splits_ : int
+        The number of cross-validation splits (folds/iterations).
+
+    refit_time_ : float
+        Seconds used for refitting the best model on the whole dataset.
+
+        This is present only if ``refit`` is not False.
+
+    multimetric_ : bool
+        Whether or not the scorers compute several metrics.
+
+    classes_ : ndarray of shape (n_classes,)
+        The classes labels. This is present only if ``refit`` is specified and
+        the underlying estimator is a classifier.
+
+    n_features_in_ : int
+        Number of features seen during :term:`fit`. Only defined if
+        `best_estimator_` is defined (see the documentation for the `refit`
+        parameter for more details) and that `best_estimator_` exposes
+        `n_features_in_` when fit.
+
+    feature_names_in_ : ndarray of shape (`n_features_in_`,)
+        Names of features seen during :term:`fit`. Only defined if
+        `best_estimator_` is defined (see the documentation for the `refit`
+        parameter for more details) and that `best_estimator_` exposes
+        `feature_names_in_` when fit.
+
+    See Also
+    --------
+    ParameterGrid : Generates all the combinations of a hyperparameter grid.
+    train_test_split : Utility function to split the data into a development
+        set usable for fitting a GridSearchCV instance and an evaluation set
+        for its final evaluation.
+    sklearn.metrics.make_scorer : Make a scorer from a performance metric or
+        loss function.
+    """
+
+    _tags = {
+        "X_inner_mtype": ["nested_univ", "numpy3D"],
+        "y_inner_mtype": ["numpy2D"],
+        "capability:multivariate": True,
+        "capability:multioutput": True,
+        "capability:unequal_length": True,
+        "capability:missing_values": True,
+        "capability:multithreading": True,
+        "capability:predict_proba": True,
+    }
+
+    def __init__(
+        self,
+        estimator,
+        param_grid,
+        scoring=None,
+        n_jobs=None,
+        refit=True,
+        cv=None,
+        verbose=0,
+        pre_dispatch="2*n_jobs",
+        error_score=np.nan,
+        return_train_score=False,
+        tune_by_variable=False,
+    ):
+        self.estimator = estimator
+        self.param_grid = param_grid
+        self.scoring = scoring
+        self.n_jobs = n_jobs
+        self.refit = refit
+        self.cv = cv
+        self.verbose = verbose
+        self.pre_dispatch = pre_dispatch
+        self.error_score = error_score
+        self.return_train_score = return_train_score
+        self.tune_by_variable = tune_by_variable
+
+        super().__init__()
+
+        DELEGATED_PARAMS = [
+            "estimator",
+            "param_grid",
+            "scoring",
+            "n_jobs",
+            "refit",
+            "cv",
+            "verbose",
+            "pre_dispatch",
+            "error_score",
+            "return_train_score",
+        ]
+
+        gcsvargs = {k: getattr(self, k) for k in DELEGATED_PARAMS}
+
+        self.estimator_ = GridSearchCV(**gcsvargs)
+
+        if self.tune_by_variable:
+            self.set_tags(**{"capability:multioutput": False})
+
+    def _fit(self, X, y):
+        """Fit time series classifier to training data.
+
+        private _fit containing the core logic, called from fit
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+            3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "pd-multiindex:":
+            pd.DataFrame with columns = variables,
+            index = pd.MultiIndex with first level = instance indices,
+            second level = time indices
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : guaranteed to be of a type in self.get_tag("y_inner_mtype")
+            1D iterable, of shape [n_instances]
+            or 2D iterable, of shape [n_instances, n_dimensions]
+            class labels for fitting
+            if self.get_tag("capaility:multioutput") = False, guaranteed to be 1D
+            if self.get_tag("capaility:multioutput") = True, guaranteed to be 2D
+
+        Returns
+        -------
+        self : Reference to self.
+        """
+        if y.shape[1] == 1:
+            y = y.flatten()
+
+        estimator = self._get_delegate()
+        estimator.fit(X=X, y=y)
+
+        fitted_param_names = [
+            "cv_results_",
+            "best_estimator_",
+            "best_score_",
+            "best_params_",
+            "best_index_",
+            "scorer_",
+            "n_splits_",
+            "refit_time_",
+            "multimetric_",
+            "classes_",
+        ]
+
+        for p in fitted_param_names:
+            if hasattr(estimator, p):
+                val = getattr(estimator, p)
+                setattr(self, p, val)
+
+        return self
+
+    def _predict(self, X):
+        """Predict labels for sequences in X.
+
+        private _predict containing the core logic, called from predict
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
+
+        Parameters
+        ----------
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+
+        Returns
+        -------
+        y : 1D np.array of int, of shape [n_instances] - predicted class labels
+            indices correspond to instance indices in X
+        """
+        estimator = self._get_delegate()
+        y_pred = estimator.predict(X=X)
+        if y_pred.ndim == 1:
+            y_pred = y_pred.reshape(-1, 1)
+        return y_pred
+
+    # the delegate is an sklearn estimator and it does not have get_fitted_params
+    # therefore we have to override _get_fitted_params from the delegator,
+    # which would otherwise call it
+    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 with str keys
+            fitted parameters, keyed by names of fitted parameter
+        """
+        return {}
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            For classifiers, a "default" set of parameters should be provided for
+            general testing, and a "results_comparison" set for comparing against
+            previously recorded results if the general set does not produce suitable
+            probabilities to compare against.
+
+        Returns
+        -------
+        params : dict or list of dict, default={}
+            Parameters to create testing instances of the class.
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`.
+        """
+        from sklearn.gaussian_process.kernels import RBF, DotProduct
+        from sklearn.metrics import accuracy_score
+
+        from sktime.classification.kernel_based import TimeSeriesSVC
+        from sktime.dists_kernels import AggrDist
+
+        mean_eucl_tskernel = AggrDist(DotProduct())
+        mean_rbf_tskernel = AggrDist(RBF())
+
+        param1 = {
+            "estimator": TimeSeriesSVC(kernel=mean_rbf_tskernel, probability=True),
+            "param_grid": {"C": [0.1, 1]},
+        }
+
+        param2 = {
+            "estimator": TimeSeriesSVC(kernel=mean_eucl_tskernel, probability=True),
+            "param_grid": {"kernel__transformer": [DotProduct(), RBF()]},
+            "scoring": accuracy_score,
+        }
+
+        return [param1, param2]
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/base/_fh.py b/sktime/forecasting/base/_fh.py
index 7a71c6e9a..dcfd619f0 100644
--- a/sktime/forecasting/base/_fh.py
+++ b/sktime/forecasting/base/_fh.py
@@ -720,15 +720,15 @@ class ForecastingHorizon:
             fh_idx = pd.Index(self.to_absolute_index(_cutoff))
 
         if cutoff is not None and isinstance(y_index, pd.MultiIndex):
-            y_inst_idx = y_index.droplevel(-1).unique()
+            y_inst_idx = y_index.droplevel(-1).unique().sort_values()
             if isinstance(y_inst_idx, pd.MultiIndex) and sort_by_time:
-                fh_list = [x + (y,) for x in y_inst_idx for y in fh_idx]
-            elif isinstance(y_inst_idx, pd.MultiIndex) and not sort_by_time:
                 fh_list = [x + (y,) for y in fh_idx for x in y_inst_idx]
+            elif isinstance(y_inst_idx, pd.MultiIndex) and not sort_by_time:
+                fh_list = [x + (y,) for x in y_inst_idx for y in fh_idx]
             elif sort_by_time:  # and not isinstance(y_inst_idx, pd.MultiIndex):
-                fh_list = [(x, y) for x in y_inst_idx for y in fh_idx]
-            else:  # not sort_by_time and not isinstance(y_inst_idx, pd.MultiIndex):
                 fh_list = [(x, y) for y in fh_idx for x in y_inst_idx]
+            else:  # not sort_by_time and not isinstance(y_inst_idx, pd.MultiIndex):
+                fh_list = [(x, y) for x in y_inst_idx for y in fh_idx]
 
             fh_idx = pd.Index(fh_list)
 
diff --git a/sktime/forecasting/compose/_reduce.py b/sktime/forecasting/compose/_reduce.py
index 460080070..2101ced93 100644
--- a/sktime/forecasting/compose/_reduce.py
+++ b/sktime/forecasting/compose/_reduce.py
@@ -184,13 +184,13 @@ def _sliding_window_transform(
         Xt = Zt[:, :, :window_length]
     # Pre-allocate array for sliding windows.
     # If the scitype is tabular regression, we have to convert X into a 2d array.
-    if scitype == "tabular-regressor":
-        if transformers is not None:
-            return yt, Xt
-        else:
-            return yt, Xt.reshape(Xt.shape[0], -1)
-    else:
-        return yt, Xt
+    if scitype == "tabular-regressor" and transformers is None:
+        Xt = Xt.reshape(Xt.shape[0], -1)
+
+    assert Xt.ndim == 2 or Xt.ndim == 3
+    assert yt.ndim == 2
+
+    return yt, Xt
 
 
 class _Reducer(_BaseWindowForecaster):
@@ -520,17 +520,16 @@ class _DirectReducer(_Reducer):
 
             if self.transformers_ is not None:
                 fh_rel = fh.to_relative(self.cutoff)
-                yt = _cut_df(yt, n_timepoints - fh_rel[i] + 1)
-                Xt = _cut_df(Xt, n_timepoints - fh_rel[i] + 1, type="head")
-                estimator.fit(Xt, yt)
+                Xt_cut = _cut_df(Xt, n_timepoints - fh_rel[i] + 1, type="head")
+                yt_cut = _cut_df(yt, n_timepoints - fh_rel[i] + 1)
+            elif self.windows_identical is True or (fh_rel[i] - 1) == 0:
+                Xt_cut = Xt
+                yt_cut = yt[:, i]
             else:
-                if self.windows_identical is True:
-                    estimator.fit(Xt, yt[:, i])
-                else:
-                    if (fh_rel[i] - 1) == 0:
-                        estimator.fit(Xt, yt[:, i])
-                    else:
-                        estimator.fit(Xt[: -(fh_rel[i] - 1)], yt[: -(fh_rel[i] - 1), i])
+                Xt_cut = Xt[: -(fh_rel[i] - 1)]
+                yt_cut = yt[: -(fh_rel[i] - 1), i]
+
+            estimator.fit(Xt_cut, yt_cut)
             self.estimators_.append(estimator)
         return self
 
diff --git a/sktime/forecasting/model_selection/_tune.py b/sktime/forecasting/model_selection/_tune.py
index 05ff5c5be..0a68e35b9 100644
--- a/sktime/forecasting/model_selection/_tune.py
+++ b/sktime/forecasting/model_selection/_tune.py
@@ -410,7 +410,7 @@ class ForecastingGridSearchCV(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"
@@ -737,7 +737,7 @@ class ForecastingRandomizedSearchCV(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"
@@ -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/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/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/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/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)