Repository: sktime/sktime

Issue:
[BUG] CV calculation is incorrect in ADICVTransformer class
**Describe the bug**
Coefficient of Variation (CV) should be calculated as the sigma of non-zero values squared over the mu of non-zero values squared; however, in the implementation here, it's being calculated as the variance of non-zero values (Which should correctly correspond to sigma squared) over the number of non-zero values, making the denominator incorrect.

This blog post contains a short description of the calculations and approach: https://frepple.com/blog/demand-classification/

Source code reference:
https://github.com/sktime/sktime/blob/cd769b1769be67888213bfff522e1ea86f71eaef/sktime/transformations/series/adi_cv.py#L167

Discord reference:
https://discord.com/channels/1075852648688930887/1075852648948973642/1261316663991861273

**To Reproduce**
Following the blog post example, here's the code to produce an incorrect ADI, CV, and Classification:

```python
import numpy as np
import pandas as pd
from datetime import datetime
from sktime.transformations.series.adi_cv import ADICVTransformer

transformer = ADICVTransformer()
periods = pd.date_range(start=datetime(2023, 1, 1), end=datetime(2023, 12, 1), freq='MS')
series = [6, 5, 0, 9, 2, 0, 14, 0, 0, 21, 0, 17]
df = pd.DataFrame({'series': series}, index=periods)

adi_cv = transformer.fit_transform(df)
print(adi_cv)
```

output:
|adi	|cv2	|class|
|:----|:-----|:----|
|2.0	|6.897959	|lumpy

**Expected behavior**
The output should look like this:
|adi	|cv2	|class|
|:----|:-----|:----|
|1.71|0.37	|intermittent


**Additional context**

**Versions**
[c:\Users\sbhobbes\Anaconda3\envs\apf\Lib\site-packages\statsforecast\core.py:27](file:///C:/Users/sbhobbes/Anaconda3/envs/apf/Lib/site-packages/statsforecast/core.py:27): TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from tqdm.autonotebook import tqdm

System:
    python: 3.11.7 | packaged by conda-forge | (main, Dec 23 2023, 14:27:59) [MSC v.1937 64 bit (AMD64)]
executable: [c:\Users\sbhobbes\Anaconda3\envs\apf\python.exe](file:///C:/Users/sbhobbes/Anaconda3/envs/apf/python.exe)
   machine: Windows-10-10.0.19042-SP0

Python dependencies:
          pip: 23.3.2
       sktime: 0.30.2
      sklearn: 1.4.0
       skbase: 0.8.1
        numpy: 1.26.3
        scipy: 1.12.0
       pandas: 2.2.2
   matplotlib: 3.8.2
       joblib: 1.3.2
        numba: 0.59.1
  statsmodels: 0.14.1
     pmdarima: 2.0.4
statsforecast: 1.7.5
      tsfresh: None
      tslearn: None
        torch: 2.3.1+cpu
   tensorflow: None
tensorflow_probability: None

</details>

<!-- Thanks for contributing! -->


Discussion:
sbhobbes: @all-contributors please add @sbhobbes for <code>
sbhobbes: @all-contributors please add @sbhobbes for :code:
sbhobbes: @all-contributors please add @sbhobbes for code
sbhobbes: @shlok191 I would appreciate any feedback/thoughts you have here as well.  Thanks!
shlok191: @sbhobbes, looks good to me. Thanks a lot again. I remember being confused about the exact ADI definition as well...There really is no 1 agreed upon definition I could find online. I think providing folks with multiple options is really helpful! :)
shlok191: Everything looks good to me too! @fkiraly
sbhobbes: @fkiraly Yeah, I'd like to make additional unit tests related to the new parameter for adi_trim_handling, I just got pulled onto a different project this week at work to help another team meet a deadline, so I haven't had a chance to get back into it. I should have something more added by the end of the weekend. 

If you want to pull it out of draft before that to get the fix into production, I'm OK with that and will do a new PR for the additional tests. 

Thank you both again!
fkiraly: the next release is later than the weekend, so no need to rush for now.
sbhobbes: > the next release is later than the weekend, so no need to rush for now - weekend is fine!

I added 2 parametrized unit tests for different expected values based on the different calculations.  I'm good with the changes, but please feel free to let me know if you or @shlok191 want me to do anything else.  Thanks again!

Patch:
diff --git a/sktime/transformations/series/adi_cv.py b/sktime/transformations/series/adi_cv.py
index c9a4f59cf..d37297048 100644
--- a/sktime/transformations/series/adi_cv.py
+++ b/sktime/transformations/series/adi_cv.py
@@ -1,6 +1,7 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Feature transformer that returns features of time series including categories."""
 
+import numpy as np
 import pandas as pd
 
 from sktime.transformations.base import BaseTransformer
@@ -15,27 +16,27 @@ class ADICVTransformer(BaseTransformer):
     ``"smooth"``, ``"erratic"``, ``"intermittent"``, or ``"lumpy"``.
 
     The labels are based on the Average Demand Interval (ADI) and the Coefficient of
-    Variation (CV2) of the time series, using simple thresholding.
-    Optionally, the transformer can return the ADI and CV2 values as well.
+    Variation squared (CV^2) of the time series, using simple thresholding.
+    Optionally, the transformer can return the ADI and CV^2 values as well.
 
     Let :math:`x_t` be the value of the time series at times :math:`t = 1, 2, \dots, T`,
     and let :math:`N` be the number of non-zero values in the time series, i.e.,
     :math:`N = \text{card}\{t: x_t \neq 0\}`.
 
-    The ADI and CV2 are calculated as follows:
+    The ADI and CV^2 are calculated as follows:
 
-    1. Average Demand Interval (ADI): The average time period between
-      time periods with non-zero demands, mathematically defined as:
+    1. Average Demand Interval (ADI): The average number of periods between
+    periods with non-zero demand, mathematically defined as:
 
-    .. math:: ADI = \frac{T}{N - 1}
+    .. math:: ADI = \frac{T}{N}
 
-    2. Variance (CV2): Variance calculated on non-zero values
+    2. Coefficient of Variation squared (CV^2): CV is calculated on non-zero values
     in the time series. Mathematical definition is:
 
     .. math:: CV2 = \frac{1}{N}\sum_{t=1}^{T} x_t^2 - \left(\frac{1}{N}\sum_{t=1}^{T} x_t\right)^2
 
     3. Class: Classification of time series on basis of ADI threshold
-    and CV2 threshold.
+    and CV^2 threshold.
 
     For thresholds ``ADI_threshold`` and ``CV2_threshold``, the classes are:
 
@@ -59,6 +60,15 @@ class ADICVTransformer(BaseTransformer):
     features : list[str] | None (default = ['adi', 'cv2', 'class'])
         Specifies all of the feature values to be calculated
 
+    adi_trim_handling: string (default = pool)
+        Specifies the method for reconciling leading/trailing zeros in the series.
+        Allowable values are ('pool', 'trim', and 'ignore'), corresponding to the
+        following treatment of leading/trailing zeros:
+        pool => L / N
+        trim => (Last N - First N) / N
+        ignore => L / (N - 1)
+            - where L is the set of All Observations, and N is the set of Non-Zero observations
+
     Examples
     --------
     >>> from sktime.transformations.series.adi_cv import ADICVTransformer
@@ -90,8 +100,14 @@ class ADICVTransformer(BaseTransformer):
         "maintainers": ["shlok191"],
     }
 
-    def __init__(self, features=None, adi_threshold=1.32, cv_threshold=0.49):
-        """Initialize the transformer and processes any provided parameters.
+    def __init__(
+        self,
+        features=None,
+        adi_threshold=1.32,
+        cv_threshold=0.49,
+        adi_trim_handling="pool",
+    ):
+        """Initialize the transformer and process any provided parameters.
 
         Parameters
         ----------
@@ -104,6 +120,10 @@ class ADICVTransformer(BaseTransformer):
         cv_threshold : float, optional
             Threshold for Variance. Defaults to 0.49.
 
+        adi_trim_handling : str, optional
+            Determines how leading and trailing zeros are handled for ADI.
+            Defaults to 'pool'
+
         Raises
         ------
             ValueError: If features is provided and does not
@@ -112,6 +132,7 @@ class ADICVTransformer(BaseTransformer):
         self.adi_threshold = adi_threshold
         self.cv_threshold = cv_threshold
         self.features = features
+        self.adi_trim_handling = adi_trim_handling
 
         self.features_internal = features
 
@@ -153,21 +174,39 @@ class ADICVTransformer(BaseTransformer):
             during initialization. Specifically, the columns include (by default):
 
                 1. Average Demand Interval (ADI)
-                2. Variance (CV2)
+                2. Coefficient of Variation squared (CV2)
                 3. categorical class
         """
-        X_non_zero = X.to_numpy().nonzero()
-        X_non_zero = X.iloc[X_non_zero]
+        X_non_zero_indices = X.to_numpy().nonzero()[0]
+        X_non_zero_values = X.iloc[X_non_zero_indices].values
 
-        # Calculating ADI value based on formula from paper
-        adi_value = len(X) / (len(X_non_zero) - 1)
+        # Calculating ADI value based on adi_trim_handling parameter
+        n_non_zero_observations = len(X_non_zero_indices)
+        n_observations = len(X)
 
-        # Calculating variance for all non-zero values
-        variance = X_non_zero.var().iloc[0]
-        cv2_value = variance / len(X_non_zero)
+        if self.adi_trim_handling == "trim":
+            first_non_zero = X_non_zero_indices[0]
+            last_non_zero = X_non_zero_indices[-1]
 
-        # Calculating the class type
+            adi_numerator = last_non_zero - first_non_zero
+            adi_denominator = n_non_zero_observations - 1
+
+        elif self.adi_trim_handling == "pool":
+            adi_numerator = n_observations
+            adi_denominator = n_non_zero_observations
+
+        elif self.adi_trim_handling == "ignore":
+            adi_numerator = n_observations
+            adi_denominator = n_non_zero_observations - 1
+
+        adi_value = adi_numerator / adi_denominator
 
+        # Calculating coefficient of variation squared for all non-zero values
+        sigma = np.std(X_non_zero_values)
+        mu = np.mean(X_non_zero_values)
+        cv2_value = (sigma / mu) ** 2
+
+        # Calculating the class type
         adi_low = adi_value <= self.adi_threshold
         cv2_low = cv2_value <= self.cv_threshold
         adi_high = not adi_low
@@ -175,13 +214,10 @@ class ADICVTransformer(BaseTransformer):
 
         if adi_low and cv2_low:
             class_type = "smooth"
-
         elif adi_low and cv2_high:
             class_type = "erratic"
-
         elif adi_high and cv2_low:
             class_type = "intermittent"
-
         elif adi_high and cv2_high:
             class_type = "lumpy"
 
@@ -232,6 +268,7 @@ class ADICVTransformer(BaseTransformer):
             "features": ["adi", "class"],
             "adi_threshold": 1.5,
             "cv_threshold": 0.2,
+            "adi_trim_handling": "trim",
         }
 
         return [params1, params2, params3]


Test patch:
diff --git a/sktime/transformations/series/tests/test_adi_cv.py b/sktime/transformations/series/tests/test_adi_cv.py
index ac08d1f3b..a0f11218a 100644
--- a/sktime/transformations/series/tests/test_adi_cv.py
+++ b/sktime/transformations/series/tests/test_adi_cv.py
@@ -9,14 +9,16 @@ from sktime.tests.test_switch import run_test_for_class
 from sktime.transformations.series.adi_cv import ADICVTransformer
 
 
-def _generate_smooth_series(size: int = 750):
+def _generate_smooth_series(size: int = 750, seed: int = 42):
     """Generates a demand time series of the "smooth" category.
-
     Parameters
     ----------
     size : int, optional
         The size of the generated time series, by default 750
 
+    seed : int, optional
+        The random seed state for the randomization algorithm.
+
     Returns
     -------
     Pandas.Series
@@ -25,19 +27,22 @@ def _generate_smooth_series(size: int = 750):
     # Generating a smooth series, we keep variance low by keeping the
     # standard deviation to 0.25. Denoted as scale according to numpy docs
 
-    smooth_series = np.random.normal(loc=10, scale=0.25, size=size)
+    np.random.seed(seed)
+    smooth_series = np.random.normal(loc=10, scale=0.25, size=12)
 
     return pd.Series(smooth_series)
 
 
-def _generate_erratic_series(size: int = 750):
+def _generate_erratic_series(size: int = 750, seed: int = 42):
     """Generates a demand time series of the "erratic" category.
-
     Parameters
     ----------
     size : int, optional
         The size of the generated time series, by default 750
 
+    seed : int, optional
+        The random seed state for the randomization algorithm.
+
     Returns
     -------
     Pandas.Series
@@ -46,19 +51,22 @@ def _generate_erratic_series(size: int = 750):
     # Generating an erratic series, we keep variance high by keeping the
     # standard deviation to 1 and then squaring the values
 
-    erratic_series = np.random.normal(loc=10, scale=2.5, size=size) ** 2
+    np.random.seed(seed)
+    erratic_series = np.random.uniform(low=1, high=100, size=size) ** 2
 
     return pd.Series(erratic_series)
 
 
-def _generate_intermittent_series(size: int = 750):
+def _generate_intermittent_series(size: int = 750, seed: int = 42):
     """Generates a demand time series of the "intermittent" category.
-
     Parameters
     ----------
     size : int, optional
         The size of the generated time series, by default 750
 
+    seed : int, optional
+        The random seed state for the randomization algorithm.
+
     Returns
     -------
     Pandas.Series
@@ -67,6 +75,7 @@ def _generate_intermittent_series(size: int = 750):
     # Generating an intermittent series, we keep ADI high by
     # setting only 10% of all values to non-zero values
 
+    np.random.seed(seed)
     intermittent_series = np.zeros(shape=(size,))
     non_zero_indices = np.random.choice(size, size=size // 10, replace=False)
 
@@ -75,14 +84,16 @@ def _generate_intermittent_series(size: int = 750):
     return pd.Series(intermittent_series)
 
 
-def _generate_lumpy_series(size: int = 750):
+def _generate_lumpy_series(size: int = 750, seed: int = 42):
     """Generates a demand time series of the "lumpy" category.
-
     Parameters
     ----------
     size : int, optional
         The size of the generated time series, by default 750
 
+    seed : int, optional
+        The random seed state for the randomization algorithm.
+
     Returns
     -------
     Pandas.Series
@@ -91,10 +102,11 @@ def _generate_lumpy_series(size: int = 750):
     # Generating a lumpy series, we keep ADI high by
     # setting only 10% of all values to non-zero values
 
+    np.random.seed(seed)
     lumpy_series = np.zeros(shape=(size,))
     non_zero_indices = np.random.choice(size, size=size // 10, replace=False)
 
-    lumpy_series[non_zero_indices] = np.random.normal(10, 2.5, size=size // 10) ** 2
+    lumpy_series[non_zero_indices] = np.random.uniform(1, 100, size=size // 10) ** 2
 
     return pd.Series(lumpy_series)
 
@@ -105,28 +117,59 @@ def _generate_lumpy_series(size: int = 750):
     reason="run test only if softdeps are present and incrementally (if requested)",
 )
 @pytest.mark.parametrize(
-    "series_generator, expected_class",
+    "demand_series, expected_adi, expected_cv, expected_class, adi_trim_handling",
     [
-        (_generate_smooth_series, "smooth"),
-        (_generate_erratic_series, "erratic"),
-        (_generate_intermittent_series, "intermittent"),
-        (_generate_lumpy_series, "lumpy"),
+        (_generate_smooth_series(), 1.0, 0.0, "smooth", "pool"),
+        (_generate_erratic_series(), 1.0, 0.82, "erratic", "pool"),
+        (_generate_intermittent_series(), 10.0, 0.0, "intermittent", "pool"),
+        (_generate_lumpy_series(), 10.0, 0.64, "lumpy", "pool"),
+        (_generate_intermittent_series(), 9.95, 0.0, "intermittent", "trim"),
+        (_generate_erratic_series(), 1.0, 0.82, "erratic", "ignore"),
+        ([10, 9, 10, 11, 12, 10, 9, 10, 8, 9, 10, 10], 1.0, 0.01, "smooth", "pool"),
+        ([10, 20, 5, 25, 20, 5, 50, 35, 30, 100, 35, 45], 1.0, 0.62, "erratic", "pool"),
+        ([6, 5, 0, 9, 2, 0, 14, 0, 0, 21, 0, 17], 1.71, 0.37, "intermittent", "pool"),
+        ([1, 0, 0, 50, 0, 0, 200, 0, 0, 0, 0, 100], 3.0, 0.70, "lumpy", "pool"),
     ],
 )
-def test_adi_cv_extractor(series_generator, expected_class):
+def test_adi_cv_extractor(
+    demand_series, expected_adi, expected_cv, expected_class, adi_trim_handling
+):
     """
     Runs a PyTest for all 4 demand time series categories.
 
     Parameters
     ----------
-    series_generator
-        A function that generates a time series in the Pandas Series format.
+    demand_series
+        A series of discrete integer values representing a demand pattern. The series
+        can be either static or randomly generated.
+
+    expected_adi
+        The expected Average Demand Interval (ADI) to be calculated by the transformer.
+
+    expected_cv
+        The expected Coefficient of Variation squared (CV^2) to be
+        calculated by the transformer.
 
     expected_class
-        The expected class to be predicted by the transformer.
+        The expected class to be returned by the transformer based on the
+        calculated ADI and CV^2 values.
     """
-    series = series_generator()
-    transformer = ADICVTransformer()
+    series = pd.Series(demand_series)
+    transformer = ADICVTransformer(adi_trim_handling=adi_trim_handling)
 
     df = transformer.fit_transform(series)
-    assert df["class"].iloc[0] == expected_class
+
+    assert np.round(df.loc[0, "adi"], 2) == expected_adi, (
+        f'The expected ADI '
+        f'of {expected_adi} does not match the actual ADI of {df.loc[0, "adi"]}.'
+    )
+
+    assert np.round(df.loc[0, "cv2"], 2) == expected_cv, (
+        f'The expected CV '
+        f'of {expected_cv} does not match the actual CV of {df.loc[0, "cv2"]}.'
+    )
+
+    assert df.loc[0, "class"] == expected_class, (
+        f'The expected class '
+        f'of {expected_class} does not match the actual class of {df.loc[0, "class"]}.'
+    )