Repository: sktime/sktime
Pull request: [BUG] Fix `ConformalIntervals` update does not update residuals_matrix
#### Reference Issues/PRs
Fixes #3419 

#### What does this implement/fix? Explain your changes.

As described in the bug, the `update` routine does not update the `residuals_matrix` which means the conformal intervals will be incorrect on the updated data. This PR computes the sliding residuals on the remaining part of the `residuals_matrix` in the case of an update.

#### Does your contribution introduce a new dependency? If yes, which one?

No

#### What should a reviewer concentrate their feedback on?

This PR adds new rows and new columns populated for the new rows, but does not extend the existing rows as it stands. Is this an acceptable approximation? The way it is used this will only matter where old rows are not extending the length of the forecasting horizon. Happy to explain this aspect in more detail if this doesn't make sense. 

#### PR checklist

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/alan-turing-institute/sktime/blob/main/.all-contributorsrc).
- [ ] Optionally, I've updated sktime's [CODEOWNERS](https://github.com/alan-turing-institute/sktime/blob/main/CODEOWNERS) to receive notifications about future changes to these files.
- [x] I've added unit tests and made sure they pass locally.
- [x] The PR title starts with either [ENH], [MNT], [DOC], or [BUG] indicating whether the PR topic is related to enhancement, maintenance, documentation, or bug.

##### For new estimators
- [ ] I've added the estimator to the online documentation.
- [ ] I've updated the existing example notebooks or provided a new one to showcase how my estimator works.



Problem statement:
[BUG] `ConformalIntervals` update does not update residuals_matrix
**Describe the bug**

`ConformalIntervals` update only updates the forecaster, not the residuals_matrix for predicting intervals or quantiles.


**To Reproduce**

`_update` is defined as:

```
    def _update(self, y, X=None, update_params=True):
        self.forecaster_.update(y, X, update_params=update_params)
        return self
```

which is only updating the forecaster and no the residuals_matrix. Since the forecaster is already fitted when `update` is called, the residuals_matrix is set to the fitted one with:

```
  if self.fh_early_:
      residuals_matrix = self.residuals_matrix_
```

**Expected behavior**

The residuals_matrix should be updated with the new data.


<!-- Thanks for contributing! -->

Hints from issue discussion:
Hm, I think both options are valid, but updating the residuals matrix should probably be the default.
Should we introduce a parameter `update_residuals` or similar?
> Hm, I think both options are valid, but updating the residuals matrix should probably be the default. Should we introduce a parameter `update_residuals` or similar?

Yes that makes sense to me (with `True` as default). Just figuring out the cleanest way to extend it without recalculating the whole thing (since there is a good chance a large portion of it has already been calculated).
I think a pain point that some forecaster fitted states have already been reached in `fit`, so there's a tradeoff between memorizing (memory use) and recomputing (compute use)

Commit message: Only update residuals matrix if update_params is true

Patch to review:
diff --git a/sktime/forecasting/conformal.py b/sktime/forecasting/conformal.py
index ce61126b5bb..27ebabb05e7 100644
--- a/sktime/forecasting/conformal.py
+++ b/sktime/forecasting/conformal.py
@@ -107,7 +107,7 @@ def __init__(
         self,
         forecaster,
         method="empirical",
-        initial_window=1,
+        initial_window=None,
         sample_frac=None,
         verbose=False,
         n_jobs=None,
@@ -162,7 +162,16 @@ def _predict(self, fh, X=None):
 
     def _update(self, y, X=None, update_params=True):
         self.forecaster_.update(y, X, update_params=update_params)
-        return self
+
+        if update_params and len(y.index.difference(self.residuals_matrix_.index)) > 2:
+            self.residuals_matrix_ = self._compute_sliding_residuals(
+                y,
+                X,
+                self.forecaster_,
+                self.initial_window,
+                self.sample_frac,
+                update=True,
+            )
 
     def _predict_interval(self, fh, X=None, coverage=None):
         """Compute/return prediction quantiles for a forecast.
@@ -289,7 +298,12 @@ def _parse_initial_window(self, y, initial_window=None):
         n_samples = len(y)
 
         if initial_window is None:
-            initial_window = max(10, 0.1 * n_samples)
+            if int(floor(0.1 * n_samples)) > 10:
+                initial_window = int(floor(0.1 * n_samples))
+            elif n_samples > 10:
+                initial_window = 10
+            else:
+                initial_window = n_samples - 1
 
         initial_window_type = np.asarray(initial_window).dtype.kind
 
@@ -317,7 +331,9 @@ def _parse_initial_window(self, y, initial_window=None):
 
         return n_initial_window
 
-    def _compute_sliding_residuals(self, y, X, forecaster, initial_window, sample_frac):
+    def _compute_sliding_residuals(
+        self, y, X, forecaster, initial_window, sample_frac, update=False
+    ):
         """Compute sliding residuals used in uncertainty estimates.
 
         Parameters
@@ -336,9 +352,12 @@ def _compute_sliding_residuals(self, y, X, forecaster, initial_window, sample_fr
             initial window.
             If None, the value is set to the larger of 0.1*len(y) and 10
         sample_frac : float
-            For speeding up computing of residuals matrix.
+            for speeding up computing of residuals matrix.
             sample value in range (0, 1) to obtain a fraction of y indices to
             compute residuals matrix for
+        update : bool
+            Whether residuals_matrix has been calculated previously and just
+            needs extending. Default = False
         Returns
         -------
         residuals_matrix : pd.DataFrame, row and column index = y.index[initial_window:]
@@ -351,11 +370,23 @@ def _compute_sliding_residuals(self, y, X, forecaster, initial_window, sample_fr
 
         n_initial_window = self._parse_initial_window(y, initial_window=initial_window)
 
-        y_index = y.iloc[n_initial_window:].index
+        full_y_index = y.iloc[n_initial_window:].index
 
-        residuals_matrix = pd.DataFrame(columns=y_index, index=y_index, dtype="float")
+        residuals_matrix = pd.DataFrame(
+            columns=full_y_index, index=full_y_index, dtype="float"
+        )
+
+        if update and hasattr(self, "residuals_matrix_"):
+            remaining_y_index = full_y_index.difference(self.residuals_matrix_.index)
+            if len(remaining_y_index) != len(full_y_index):
+                residuals_matrix.loc[
+                    self.residuals_matrix_.index, self.residuals_matrix_.columns
+                ] = self.residuals_matrix_
+            y_index = remaining_y_index
+        else:
+            y_index = full_y_index
 
-        if sample_frac:
+        if sample_frac and len(y_index.to_series().sample(frac=sample_frac)) > 2:
             y_index = y_index.to_series().sample(frac=sample_frac)
 
         def _get_residuals_matrix_row(forecaster, y, X, id):
@@ -364,7 +395,6 @@ def _get_residuals_matrix_row(forecaster, y, X, id):
 
             X_train = get_slice(X, start=None, end=id)
             X_test = get_slice(X, start=id, end=None)
-
             forecaster.fit(y_train, X=X_train, fh=y_test.index)
 
             try:
@@ -380,7 +410,6 @@ def _get_residuals_matrix_row(forecaster, y, X, id):
             delayed(_get_residuals_matrix_row)(forecaster.clone(), y, X, id)
             for id in y_index
         )
-
         for idx, id in enumerate(y_index):
             residuals_matrix.loc[id] = all_residuals[idx]
 


Code review comments:
diff_hunk: @@ -351,11 +370,23 @@ def _compute_sliding_residuals(self, y, X, forecaster, initial_window, sample_fr
 
         n_initial_window = self._parse_initial_window(y, initial_window=initial_window)
 
-        y_index = y.iloc[n_initial_window:].index
+        full_y_index = y.iloc[n_initial_window:].index
 
-        residuals_matrix = pd.DataFrame(columns=y_index, index=y_index, dtype="float")
+        residuals_matrix = pd.DataFrame(
+            columns=full_y_index, index=full_y_index, dtype="float"
+        )
+
+        if update and hasattr(self, "residuals_matrix_"):
+            remaining_y_index = full_y_index.difference(self.residuals_matrix_.index)
+            if len(remaining_y_index) != len(full_y_index):
+                residuals_matrix.loc[
+                    self.residuals_matrix_.index, self.residuals_matrix_.columns
+                ] = self.residuals_matrix_
+            y_index = remaining_y_index
+        else:
+            y_index = full_y_index
 
-        if sample_frac:
+        if sample_frac and len(y_index.to_series().sample(frac=sample_frac)) > 2:
path: sktime/forecasting/conformal.py
text: @user1:
you are sampling twice, this can lead to the sample in the "if" satisfying the condition and the other sample not. So, should we not sample and then check the lenght of *that* sample?

@author:
Oh yes - fixed now