Feat/probabilistic ensemble

CHANGELOG.md

@@ -5,6 +5,9 @@ We do our best to avoid the introduction of breaking changes,
 but cannot always guarantee backwards compatibility. Changes that may **break code which uses a previous release of Darts** are marked with a "🔴".
 
 ## [Unreleased](https://github.com/unit8co/darts/tree/master)
+- Improvements to `EnsembleModel`:
+  - `RegressionEnsembleModel` and `NaiveEnsembleModel` can generate probabilistic forecasts, probabilistics `forecasting_models` can be sampled to train the `regression_model`, updated the documentation (stacking technique). [#1692](https://github.com/unit8co/darts/pull/#1692) by [Antoine Madrona](https://github.com/madtoinou).
+
 [Full Changelog](https://github.com/unit8co/darts/compare/0.24.0...master)
 
 ### For users of the library:

darts/models/forecasting/baselines.py

@@ -170,7 +170,7 @@ def __init__(
         Naive implementation of `EnsembleModel`
         Returns the average of all predictions of the constituent models
         """
-        super().__init__(models)
+        super().__init__(models, train_num_samples=None, train_samples_reduction=None)
 
     def fit(
         self,
@@ -200,11 +200,13 @@ def ensemble(
         self,
         predictions: Union[TimeSeries, Sequence[TimeSeries]],
         series: Optional[Sequence[TimeSeries]] = None,
+        num_samples: int = 1,
     ) -> Union[TimeSeries, Sequence[TimeSeries]]:
         def take_average(prediction: TimeSeries) -> TimeSeries:
-            series = prediction.pd_dataframe(copy=False).sum(axis=1) / len(self.models)
-            series.name = prediction.components[0]
-            return TimeSeries.from_series(series)
+            # average across the components, keep n_samples, rename components
+            return prediction.mean(axis=1).with_columns_renamed(
+                "components_mean", prediction.components[0]
+            )
 
         if isinstance(predictions, Sequence):
             return [take_average(p) for p in predictions]

darts/models/forecasting/ensemble_model.py

@@ -12,6 +12,7 @@
     LocalForecastingModel,
 )
 from darts.timeseries import TimeSeries
+from darts.utils.utils import series2seq
 
 logger = get_logger(__name__)
 
@@ -26,10 +27,24 @@ class EnsembleModel(GlobalForecastingModel):
     ----------
     models
         List of forecasting models whose predictions to ensemble
+
+        .. note::
+                if all the models are probabilistic, the `EnsembleModel` will also be probabilistic.
+        ..
+    train_num_samples
+        Number of prediction samples from each forecasting model for multi-level ensembles. The n_samples
+        dimension will be reduced using the `train_samples_reduction` method.
+    train_samples_reduction
+        If `forecasting models` are probabilistic and `train_num_samples` > 1, method used to
+        reduce the samples dimension to 1. Possible values: "mean", "median" or float value corresponding
+        to the desired quantile.
     """
 
     def __init__(
-        self, models: Union[List[LocalForecastingModel], List[GlobalForecastingModel]]
+        self,
+        models: Union[List[LocalForecastingModel], List[GlobalForecastingModel]],
+        train_num_samples: int,
+        train_samples_reduction: Union[str, float],
     ):
         raise_if_not(
             isinstance(models, list) and models,
@@ -59,7 +74,10 @@ def __init__(
 
         super().__init__()
         self.models = models
+        self.train_num_samples = train_num_samples
+        self.train_samples_reduction = train_samples_reduction
 
+    @abstractmethod
     def fit(
         self,
         series: Union[TimeSeries, Sequence[TimeSeries]],
@@ -160,20 +178,31 @@ def predict(
             verbose=verbose,
         )
 
+        # for multi-level models, forecasting models can generate arbitrary number of samples
+        if self.train_samples_reduction is None:
+            pred_num_samples = num_samples
+        else:
+            pred_num_samples = self.train_num_samples
+
         predictions = self._make_multiple_predictions(
             n=n,
             series=series,
             past_covariates=past_covariates,
             future_covariates=future_covariates,
-            num_samples=num_samples,
+            num_samples=pred_num_samples,
         )
-        return self.ensemble(predictions, series=series)
+
+        if self.train_samples_reduction is not None and self.train_num_samples > 1:
+            predictions = self._predictions_reduction(predictions)
+
+        return self.ensemble(predictions, series=series, num_samples=num_samples)
 
     @abstractmethod
     def ensemble(
         self,
         predictions: Union[TimeSeries, Sequence[TimeSeries]],
         series: Optional[Sequence[TimeSeries]] = None,
+        num_samples: int = 1,
     ) -> Union[TimeSeries, Sequence[TimeSeries]]:
         """
         Defines how to ensemble the individual models' predictions to produce a single prediction.
@@ -193,6 +222,20 @@ def ensemble(
         """
         pass
 
+    def _predictions_reduction(self, predictions: TimeSeries) -> TimeSeries:
+        """Reduce the sample dimension of the forecasting models predictions"""
+        is_single_series = isinstance(predictions, TimeSeries)
+        predictions = series2seq(predictions)
+        if self.train_samples_reduction == "median":
+            predictions = [pred.median(axis=2) for pred in predictions]
+        elif self.train_samples_reduction == "mean":
+            predictions = [pred.mean(axis=2) for pred in predictions]
+        else:
+            predictions = [
+                pred.quantile(self.train_samples_reduction) for pred in predictions
+            ]
+        return predictions[0] if is_single_series else predictions
+
     @property
     def min_train_series_length(self) -> int:
         return max(model.min_train_series_length for model in self.models)
@@ -227,5 +270,8 @@ def find_max_lag_or_none(lag_id, aggregator) -> Optional[int]:
             find_max_lag_or_none(i, agg) for i, agg in enumerate(lag_aggregators)
         )
 
-    def _is_probabilistic(self) -> bool:
+    def _models_are_probabilistic(self) -> bool:
         return all([model._is_probabilistic() for model in self.models])
+
+    def _is_probabilistic(self) -> bool:
+        return self._models_are_probabilistic()

darts/models/forecasting/regression_ensemble_model.py

@@ -19,6 +19,9 @@
 
 logger = get_logger(__name__)
 
+# arbitrary threshold to raise warning for probabilistic forecasting models
+SAMPLES_WARNING_THRESHOLD = 1e6
+
 
 class RegressionEnsembleModel(EnsembleModel):
     def __init__(
@@ -28,9 +31,11 @@ def __init__(
         ],
         regression_train_n_points: int,
         regression_model=None,
+        regression_train_num_samples: Optional[int] = 1,
+        regression_train_samples_reduction: Optional[Union[str, float]] = "median",
     ):
         """
-        Use a regression model for ensembling individual models' predictions.
+        Use a regression model for ensembling individual models' predictions using the stacking technique [1]_.
 
         The provided regression model must implement ``fit()`` and ``predict()`` methods
         (e.g. scikit-learn regression models). Note that here the regression model is used to learn how to
@@ -47,8 +52,26 @@ def __init__(
         regression_model
             Any regression model with ``predict()`` and ``fit()`` methods (e.g. from scikit-learn)
             Default: ``darts.model.LinearRegressionModel(fit_intercept=False)``
+
+            .. note::
+                if `regression_model` is probabilistic, the `RegressionEnsembleModel` will also be probabilistic.
+            ..
+        regression_train_num_samples
+            Number of prediction samples from each forecasting model to train the regression model (samples are
+            averaged). Should be set to 1 for deterministic models. Default: 1.
+        regression_train_samples_reduction
+            If `forecasting models` are probabilistic and `regression_train_num_samples` > 1, method used to
+            reduce the samples before passing them to the regression model. Possible values: "mean", "median"
+            or float value corresponding to the desired quantile. Default: "median"
+        References
+        ----------
+        .. [1] D. H. Wolpert, “Stacked generalization”, Neural Networks, vol. 5, no. 2, pp. 241–259, Jan. 1992
         """
-        super().__init__(forecasting_models)
+        super().__init__(
+            forecasting_models,
+            train_num_samples=regression_train_num_samples,
+            train_samples_reduction=regression_train_samples_reduction,
+        )
         if regression_model is None:
             regression_model = LinearRegressionModel(
                 lags=None, lags_future_covariates=[0], fit_intercept=False
@@ -69,6 +92,52 @@ def __init__(
             f"{regression_model.lags}",
         )
 
+        raise_if(
+            regression_train_num_samples > 1 and not self._models_are_probabilistic(),
+            "`regression_train_num_samples` is greater than 1 but the `RegressionEnsembleModel` "
+            "contains at least one non-probabilistic forecasting model.",
+            logger,
+        )
+
+        # check the reduction method
+        supported_reduction = ["mean", "median"]
+        if isinstance(regression_train_samples_reduction, float):
+            # this is already checked by `ts.quantile()`, maybe too redundant
+            raise_if(
+                regression_train_samples_reduction > 1.0
+                or regression_train_samples_reduction < 0,
+                f"`regression_train_samples_reduction` should be comprised between "
+                f"0 and 1 ({regression_train_samples_reduction}).",
+                logger,
+            )
+        elif isinstance(regression_train_samples_reduction, str):
+            raise_if(
+                regression_train_samples_reduction not in supported_reduction,
+                f"`regression_train_samples_reduction` should be one of {supported_reduction}, "
+                f"received ({regression_train_samples_reduction})",
+                logger,
+            )
+        else:
+            logger.exception(
+                f"`regression_train_samples_reduction` type not supported "
+                f"({regression_train_samples_reduction}). Must be `float` "
+                f" or one of {supported_reduction}."
+            )
+
+        if (
+            regression_train_num_samples
+            * regression_train_n_points
+            * len(forecasting_models)
+            > SAMPLES_WARNING_THRESHOLD
+        ):
+            logger.warning(
+                f"Considering the number of models present in this ensemble ({len(forecasting_models)}), "
+                f"`regression_train_n_points` ({regression_train_n_points}) and `regression_train_num_samples` "
+                f"({regression_train_num_samples}) the number of sampled values to train the regression model "
+                f"will be very large ({regression_train_num_samples*regression_train_n_points*len(forecasting_models)}"
+                f">{SAMPLES_WARNING_THRESHOLD})."
+            )
+
         self.regression_model = regression_model
         self.train_n_points = regression_train_n_points
 
@@ -101,7 +170,7 @@ def fit(
 
         raise_if(
             train_n_points_too_big,
-            "regression_train_n_points parameter too big (must be smaller or "
+            "`regression_train_n_points` parameter too big (must be smaller or "
             "equal to the number of points in training_series)",
             logger,
         )
@@ -126,9 +195,13 @@ def fit(
             series=forecast_training,
             past_covariates=past_covariates,
             future_covariates=future_covariates,
-            num_samples=1,
+            num_samples=self.train_num_samples,
         )
 
+        # component-wise reduction of the probabilistic forecasting models predictions
+        if predictions[0].n_samples > 1:
+            predictions = self._predictions_reduction(predictions)
+
         # train the regression model on the individual models' predictions
         self.regression_model.fit(
             series=regression_target, future_covariates=predictions
@@ -152,6 +225,7 @@ def ensemble(
         self,
         predictions: Union[TimeSeries, Sequence[TimeSeries]],
         series: Optional[Sequence[TimeSeries]] = None,
+        num_samples: int = 1,
     ) -> Union[TimeSeries, Sequence[TimeSeries]]:
 
         is_single_series = isinstance(series, TimeSeries) or series is None
@@ -160,7 +234,10 @@ def ensemble(
 
         ensembled = [
             self.regression_model.predict(
-                n=len(prediction), series=serie, future_covariates=prediction
+                n=len(prediction),
+                series=serie,
+                future_covariates=prediction,
+                num_samples=num_samples,
             )
             for serie, prediction in zip(series, predictions)
         ]
@@ -179,3 +256,10 @@ def extreme_lags(
     ]:
         extreme_lags_ = super().extreme_lags
         return (extreme_lags_[0] - self.train_n_points,) + extreme_lags_[1:]
+
+    def _is_probabilistic(self) -> bool:
+        """
+        A RegressionEnsembleModel is probabilistic if its regression
+        model is probabilistic (ensembling layer)
+        """
+        return self.regression_model._is_probabilistic()