Repository: sktime/sktime
Pull request: [BUG] fix `BATS` and `TBATS` `_predict_interval` interface
Fixes https://github.com/sktime/sktime/issues/4491

The `_predict_interval` interface in `_TbatsAdapter` was a bit of a mess.

There was also really strange code that triggered specifically whenever the `fh` was longer than the training data, I was unable to understand why that condition should be the same as "contains in-sample prediction intervals".

This is the code:
```python
        # If pred_int contains in-sample prediction intervals
        if len(pred_int) > len(self._y):
            len_out = len(pred_int) - len(self._y)
            # Workaround for slicing with negative index
            pred_int["idx"] = [x for x in range(-len(self._y), len_out)]
```

It explains @ngupta23's issues in #4491, but I don't get why this was there in the first place.

Either way, I just completely rewrote the interface to prediction intervals from scratch.
The logic was simple, it is producing `nan` for insample predictions.

Contains the failure case from https://github.com/sktime/sktime/issues/4491 as a test.

Problem statement:
[BUG] TBATS can not produce a forecast with intervals when FH > History Period
**Describe the bug**
<!--
A clear and concise description of what the bug is.
-->

**To Reproduce**
```python
import numpy as np
import pandas as pd
import logging
from sktime.forecasting.tbats import TBATS


# generate train set values
np.random.seed(42)
LEN_HISTORY = 50
train = pd.Series(data=np.random.randint(1, 100, LEN_HISTORY))

# train model
estimator = TBATS(
    use_box_cox=False,
    use_trend=True,
    use_damped_trend=False,
    sp=10,
    use_arma_errors=False,
)
estimator.fit(train)

# FH = Historical length works!
try:
    fcst1 = estimator.predict_interval(
        coverage=0.8, fh=np.array(range(1, LEN_HISTORY + 1))
    )
except Exception as e:
    print("FH <= History")
    logging.exception(e)

# FH > Historical length DOES NOT work
try:
    fcst2 = estimator.predict_interval(
        coverage=0.8, fh=np.array(range(1, LEN_HISTORY + 2))
    )
except Exception as e:
    print("FH > History")
    logging.exception(e)
```

```python-traceback
FH > History
ERROR:root:Length mismatch: Expected axis has 1 elements, new values have 51 elements
Traceback (most recent call last):
  File "analysis/tbats/TBATS_prediction_interval_issue.py", line 34, in <module>
    coverage=0.8, fh=np.array(range(1, LEN_HISTORY + 2))
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/sktime/forecasting/base/_base.py", line 636, in predict_interval
    pred_int = self._predict_interval(fh=fh, X=X_inner, coverage=coverage)
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/sktime/forecasting/base/adapters/_tbats.py", line 293, in _predict_interval
    _, tbats_pred_int = self._tbats_forecast_with_interval(fh, c)
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/sktime/forecasting/base/adapters/_tbats.py", line 217, in _tbats_forecast_with_interval
    pred_int = self._get_pred_int(lower=lower, upper=upper)
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/sktime/forecasting/base/adapters/_tbats.py", line 345, in _get_pred_int
    pred_int.index = fh_out.to_absolute(self.cutoff)
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/pandas/core/generic.py", line 5154, in __setattr__
    return object.__setattr__(self, name, value)
  File "pandas/_libs/properties.pyx", line 66, in pandas._libs.properties.AxisProperty.__set__
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/pandas/core/generic.py", line 564, in _set_axis
    self._mgr.set_axis(axis, labels)
  File "/home/<hidden>/anaconda3/envs/<hidden>/lib/<hidden>/site-packages/pandas/core/internals/managers.py", line 227, in set_axis
    f"Length mismatch: Expected axis has {old_len} elements, new "
ValueError: Length mismatch: Expected axis has 1 elements, new values have 51 elements
```

**Expected behavior**
Requesting Prediction Intervals should work when FH > len(history) since it works when just requesting predictions without intervals (so it there should not be a historical data insufficiency issue).


**Versions**
<details>
sktime 0.16.0 and sktime 0.17.0
</details>

<!-- Thanks for contributing! -->

Hints from issue discussion:
FYI @yarnabrina,this looks similar to code you are currently refactoring for `statsmodels`.
This seems to come from `fh`, have you encountered this before?
update: bug confirmed, python 3.10, current main

This is also not coming from `fh`, that is fine.
It is coming from preceding code, which produces `pred_int` of wrong length, as soon as the prediction index exceeds length of the previous series.

Commit message: add test case from 4491

Patch to review:
diff --git a/sktime/forecasting/base/adapters/_tbats.py b/sktime/forecasting/base/adapters/_tbats.py
index b4e8db6b945..c59c41765b9 100644
--- a/sktime/forecasting/base/adapters/_tbats.py
+++ b/sktime/forecasting/base/adapters/_tbats.py
@@ -3,7 +3,7 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Implements adapter for using tbats forecasters in sktime framework."""
 
-__author__ = ["mloning", "aiwalter", "k1m190r"]
+__author__ = ["mloning", "aiwalter", "k1m190r", "fkiraly"]
 __all__ = ["_TbatsAdapter"]
 
 import numpy as np
@@ -178,7 +178,7 @@ def _get_y_pred(self, y_in_sample, y_out_sample):
         return y_pred
 
     def _tbats_forecast(self, fh):
-        """TBATS forecast without confidence interval.
+        """TBATS point forecast adapter function.
 
         Parameters
         ----------
@@ -208,8 +208,8 @@ def _tbats_forecast(self, fh):
 
         return y_pred
 
-    def _tbats_forecast_with_interval(self, fh, conf_lev):
-        """TBATS forecast with confidence interval.
+    def _tbats_forecast_interval(self, fh, conf_lev):
+        """TBATS prediction interval forecast adapter function.
 
         Parameters
         ----------
@@ -226,42 +226,22 @@ def _tbats_forecast_with_interval(self, fh, conf_lev):
             Prediction intervals
         """
         fh = fh.to_relative(cutoff=self.cutoff)
-        len_fh = len(fh)
+        fh_out = fh.to_out_of_sample(cutoff=self.cutoff)
+        steps = fh_out.to_pandas().max()
 
-        if not fh.is_all_in_sample(cutoff=self.cutoff):
-            fh_out = fh.to_out_of_sample(cutoff=self.cutoff)
-            steps = fh_out.to_pandas().max()
-            _, tbats_ci = self._forecaster.forecast(
-                steps=steps, confidence_level=conf_lev
-            )
-            out = pd.DataFrame(tbats_ci)
-            y_out = out["mean"]  # aka tbats y_hat out of sample
-
-            # pred_int
-            lower = pd.Series(out["lower_bound"])
-            upper = pd.Series(out["upper_bound"])
-            pred_int = self._get_pred_int(lower=lower, upper=upper)
-
-            if len(fh) != len(fh_out):
-                epred_int = pd.DataFrame({"lower": nans(len_fh), "upper": nans(len_fh)})
-                epred_int.index = fh.to_absolute(self.cutoff).to_pandas()
-
-                in_pred_int = epred_int.index.isin(pred_int.index)
-                epred_int[in_pred_int] = pred_int
-                pred_int = epred_int
+        _, tbats_ci = self._forecaster.forecast(steps=steps, confidence_level=conf_lev)
+        out = pd.DataFrame(tbats_ci)
 
-        else:
-            y_out = nans(len_fh)
-            pred_int = pd.DataFrame({"lower": nans(len_fh), "upper": nans(len_fh)})
-            pred_int.index = fh.to_absolute(self.cutoff).to_pandas()
-
-        # y_pred
-        y_in_sample = pd.Series(self._forecaster.y_hat)
-        y_out_sample = pd.Series(y_out)
-        y_pred = self._get_y_pred(y_in_sample=y_in_sample, y_out_sample=y_out_sample)
-        y_pred.name = self._yname
+        # pred_int
+        lower = pd.Series(out["lower_bound"])
+        upper = pd.Series(out["upper_bound"])
+        pred_int_oos = pd.DataFrame({"lower": lower, "upper": upper})
+        pred_int_oos = pred_int_oos.iloc[fh_out.to_indexer()]
+        pred_int_oos.index = fh_out.to_absolute(self.cutoff).to_pandas()
+        full_ix = fh.to_absolute(self.cutoff).to_pandas()
+        pred_int = pred_int_oos.reindex(full_ix)
 
-        return y_pred, pred_int
+        return pred_int
 
     def _predict_interval(self, fh, X, coverage):
         """Compute/return prediction quantiles for a forecast.
@@ -317,8 +297,8 @@ def _predict_interval(self, fh, X, coverage):
                 pred_int[("Coverage", 0, "upper")] = pred_int[("Coverage", 0, "lower")]
                 continue
 
-            # tbats with CI intervals
-            _, tbats_pred_int = self._tbats_forecast_with_interval(fh, c)
+            # tbats prediction intervals
+            tbats_pred_int = self._tbats_forecast_interval(fh, c)
 
             pred_int[("Coverage", c, "lower")] = tbats_pred_int["lower"]
             pred_int[("Coverage", c, "upper")] = tbats_pred_int["upper"]
@@ -341,39 +321,6 @@ def _get_fitted_param_names(self):
         """Get names of fitted parameters."""
         return self._fitted_param_names
 
-    def _get_pred_int(self, lower, upper):
-        """Combine lower/upper bounds of pred.intervals, slice on fh.
-
-        Parameters
-        ----------
-        lower : pd.Series
-            Lower bound (can contain also in-sample bound)
-        upper : pd.Series
-            Upper bound (can contain also in-sample bound)
-
-        Returns
-        -------
-        pd.DataFrame
-            pred_int, prediction intervals (out-sample, sliced by fh)
-        """
-        pred_int = pd.DataFrame({"lower": lower, "upper": upper})
-        # Out-sample fh
-        fh_out = self.fh.to_out_of_sample(cutoff=self.cutoff)
-        # If pred_int contains in-sample prediction intervals
-        if len(pred_int) > len(self._y):
-            len_out = len(pred_int) - len(self._y)
-            # Workaround for slicing with negative index
-            pred_int["idx"] = [x for x in range(-len(self._y), len_out)]
-        # If pred_int does not contain in-sample prediction intervals
-        else:
-            pred_int["idx"] = [x for x in range(len(pred_int))]
-        pred_int = pred_int.loc[
-            pred_int["idx"].isin(fh_out.to_indexer(self.cutoff).values)
-        ]
-        pred_int.index = fh_out.to_absolute(self.cutoff).to_pandas()
-        pred_int = pred_int.drop(columns=["idx"])
-        return pred_int
-
 
 def nans(length):
     """Return a vector of NaNs, of length `length`."""
diff --git a/sktime/forecasting/tests/test_tbats.py b/sktime/forecasting/tests/test_tbats.py
new file mode 100644
index 00000000000..a60160731ed
--- /dev/null
+++ b/sktime/forecasting/tests/test_tbats.py
@@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+"""Tests for TBATS."""
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+
+__author__ = ["fkiraly", "ngupta23"]
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from sktime.forecasting.tbats import TBATS
+from sktime.utils.validation._dependencies import _check_estimator_deps
+
+
+@pytest.mark.skipif(
+    not _check_estimator_deps(TBATS, severity="none"),
+    reason="skip test if required soft dependency not available",
+)
+def test_tbats_long_fh():
+    """Test TBATS with long fh, checks for failure condition in bug #4491."""
+    np.random.seed(42)
+    LEN_HISTORY = 50
+    train = pd.Series(data=np.random.randint(1, 100, LEN_HISTORY))
+
+    # train model
+    estimator = TBATS(
+        use_box_cox=False,
+        use_trend=True,
+        use_damped_trend=False,
+        sp=10,
+        use_arma_errors=False,
+    )
+    estimator.fit(train)
+
+    # failure condition is fh being longer than training data
+    long_fh = np.array(range(1, LEN_HISTORY + 2))
+
+    fcst = estimator.predict_interval(coverage=0.8, fh=long_fh)
+    assert len(fcst) == len(long_fh)
diff --git a/sktime/forecasting/tests/test_theta.py b/sktime/forecasting/tests/test_theta.py
index f2b2fc14ae7..724f65601e3 100644
--- a/sktime/forecasting/tests/test_theta.py
+++ b/sktime/forecasting/tests/test_theta.py
@@ -1,8 +1,6 @@
 # -*- coding: utf-8 -*-
-"""Tests for ThetaForecaster.
-
+"""Tests for ThetaForecaster."""
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
-"""
 
 __author__ = ["big-o", "kejsitake"]
 


Code review comments:
diff_hunk: @@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+"""Tests for TBATS."""
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+
+__author__ = ["fkiraly", "ngupta23"]
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from sktime.forecasting.tbats import TBATS
+from sktime.utils.validation._dependencies import _check_estimator_deps
+
+
+@pytest.mark.skipif(
+    not _check_estimator_deps(TBATS, severity="none"),
+    reason="skip test if required soft dependency not available",
+)
+def test_tbats_long_fh():
+    """Test TBATS with long fh, checks for failure condition in bug #4491."""
+    np.random.seed(42)
path: sktime/forecasting/tests/test_tbats.py
text: @user1:
Do unit tests inputs need to be reproducible? If not, removing this will ensure different sample in different runs.

@author:
hm, I suppose not, unless you try to reproduce hard-coded expected outcomes.

I just copy-pasted this from @user2's example in #4491, and I wouldn't see a problem with having this line, or not having it.

Do you have a preference?

@user2:
My 2 cents on this - I feel it is better to be reproducible in tests so you can quickly distinguish actual failures rather than one that is due to the data changing from test to test. If we need to test with multiple datasets (with different conditions), it is better to parameterize the data used for testing (but still be explicit with the dataset definition).

@author:
yes, I suppose that's the main advantage.
The disadvantage being, if you have a sporadic failure - like the one with quantiles crossing on occasion - you'll never spot it.

(not arguing for either option, still on the fence)

@user1:
I think unit tests should test with as general example as possible. I like `hypothesis` a lot specifically for this, if my function takes a numeric input, I should be able to test for any (does not necessarily mean all to save resources) number, not just 5 or 6 values I happen to choose at time of writing, as almost never ever those get modified afterwards and test for other values. Apart for any general dataset, I'd also want to test special cases (for example 0 if division is involved), and for that I'd prefer hardcoded inputs, but only because they are special.

As for this example, we need to test if `tbats` support to predict intervals or not, not if it can support for the specific data points that `seed=42` would generate. Our goal is to ensure that `tbats` works, not to pass tests and get a green tick in GHA. So, ideally I would prefer to test `tbats` for any data, and hopefully quite a good amount of variation. But obviously that would take some good amount of time, so I'd agree to restrict to very restricted number of examples, even 1. But I think we should at least test estimator's generalisation for other datasets over multiple runs in different PR's.

Had you used some sort of "gold/standard dataset", to me it'd make sense to fix the inputs. Since it is not, I'd prefer to not fix it.

(I sincerely apologise if this sounds too harsh or opinionated. I am just trying to share the point of view from which I suggested. I understand there are other tests in `sktime` which also use `seed`, so it should be okay if this line is kept as well.)

@author:
Hm, good points, @user1, and thanks for the learning moment on `hypothesis` - really neat package!

I would now agree with what you say, I think "covering generalisation" is better than having the test reproducible (especially since we are not reproducing against expected results).

A counterargument to @user2 would be to

> actual failures rather than one that is due to the data changing from test to test

that is, if the data changes from test to test and that's why it fails - that *is* an actual failure, i.e., something we would want to know about!
We now fixed the funny condition of `fh` longer than the training set, but what if there are other funny conditions, such as, I don't know, `fh` length is a multiple of 3 or 7, and contains at least two negative entries.

Another way to think about this would be: the very reason why we did not spot #4491 before was the lack of variation in the test examples. We wouldn't have had the bug "hiding" in the first place if we would have had more variation.

(it also works in favour of @user2 if we are more likely to spot obscure edge failure cases than not)

@author:
(removed the random seed, following @user1's argument)

@user2:
Those are good arguments. 

Quick question: If a test were to fail with `hypothesis` on GitHub CI, would you be able to reproduce it locally for finding RCA and fixing it? If so, how would this be accomplished?


@author:
> would you be able to reproduce it locally for finding RCA and fixing it? If so, how would this be accomplished?

Excellent question! @user1, I'd also be interested in the answer!

You can of course always brute-force it and just sample 100 or 1000 times and break on exception, but probably there's a more elegant way than this.

@author:
found the relevant paragraph in the docs of `hypothesis`:
https://hypothesis.readthedocs.io/en/latest/reproducing.html

> When a test fails unexpectedly, usually due to a health check failure, Hypothesis will print out a seed that led to that failure

@user1:
> found the relevant paragraph in the docs of `hypothesis`: https://hypothesis.readthedocs.io/en/latest/reproducing.html
> 
> > When a test fails unexpectedly, usually due to a health check failure, Hypothesis will print out a seed that led to that failure

@author While this is an option, you still have to run all the test(s) with that seed for lot of iteration for it to fail once with a single case. If there is a actual system failure or something, that would help, but I interpreted @user2's question as to find the particular choice of example that led to failure. There's a simpler option to do that, as all failing tests will show the corresponding inputs. Here's an example from [Quick start guide](https://hypothesis.readthedocs.io/en/latest/quickstart.html):

> Falsifying example: test_decode_inverts_encode(s='')

Also, instead of just returning all or first failing examples, by default `hypothesis` finds the "minimal" (subject to defined [shrinking](https://hypothesis.readthedocs.io/en/latest/data.html#shrinking) strategy) example.

<details>

<summary>
Few `hypothesis` examples (nothing to do with `sktime`, just for demonstration)
</summary>

```python3
import collections.abc

import hypothesis
import hypothesis.strategies


def sample_function_1(input_list: list[int]) -> float:
    """Should fail if ``input_list`` is empty, or contains `0`."""
    input_length = len(input_list)
    minimum_absolute_value = min(abs(value) for value in input_list)

    return input_length / minimum_absolute_value


def sample_function_2(input_list: list[int]) -> float:
    """Should fail if contains `0`."""
    # handle empty case
    if not input_list:
        return 0.0

    input_length = len(input_list)
    minimum_absolute_value = min(abs(value) for value in input_list)

    return input_length / minimum_absolute_value


def sample_function_3(input_list: list[int]) -> float:
    """Should not fail."""
    # handle empty case
    if not input_list:
        return 0.0

    # handle zero case
    if 0 in input_list:
        return float("inf")

    input_length = len(input_list)
    minimum_absolute_value = min(abs(value) for value in input_list)

    return input_length / minimum_absolute_value


@user3.given(
    sample_list=hypothesis.strategies.lists(hypothesis.strategies.integers(), min_size=1)
)
def test_sample_function_1(sample_list: list[int]) -> None:
    """Useless test that always pass, unless ``sample_function_1`` fails."""
    hypothesis.assume(0 not in sample_list)

    return_value = sample_function_1(sample_list)

    assert isinstance(return_value, float)


@user3.given(sample_list=hypothesis.strategies.lists(hypothesis.strategies.integers()))
def test_sample_function_2(sample_list: list[int]) -> None:
    """Useless test that always pass, unless ``sample_function_2`` fails."""
    hypothesis.assume(0 not in sample_list)

    return_value = sample_function_2(sample_list)

    assert isinstance(return_value, float)


@user3.given(sample_list=hypothesis.strategies.lists(hypothesis.strategies.integers()))
def test_sample_function_3(sample_list: list[int]) -> None:
    """Useless test that always pass, unless ``sample_function_3`` fails."""
    return_value = sample_function_3(sample_list)

    assert isinstance(return_value, float)


@user3.given(
    sample_list=hypothesis.strategies.lists(hypothesis.strategies.integers()),
    sample_function=hypothesis.strategies.sampled_from(
        [sample_function_1, sample_function_2, sample_function_3]
    ),
)
def test_sample_function(
    sample_list: list[int], sample_function: collections.abc.Callable[[list[int]], float]
) -> None:
    """Useless test that always pass, unless ``sample_function`` fails."""
    return_value = sample_function(sample_list)

    assert isinstance(return_value, float)
```

</details>

If you run this with `pytest test_hypothesis.py` in an environment with `hypothesis[pytest]` installed, first 3 should pass and last one should fail, and it'll should show that `sample_function_1` failed with `[]`, and from that you can start debugging.

```console
E       Falsifying example: test_sample_function(
E           sample_list=[],
E           sample_function=sample_function_1,
E       )
```

@author:
nice! `hypothesis` would be a very interesting thing to show in the community colalb sessions (if you have time?), we could think about rebasing parts of the test framework on it!

Theyare also plugin extensible, so we could think about how to make the fixture generation logic (which framework-wise sits in `skbase`) compatible with `hypothesis`.

Going slightly off-topic here. Might be worth a design issue!