Repository: Ouranosinc/xclim

Issue:
Initial values for distribution fit
* xclim version: 0.19

### Description
<!--Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.-->

The scipy.stats fit function using Maximum Likelihood (`generic.fit`)  converges to poor values for some distributions (weibull_min, expon, genextreme). 

In R's `gev.fit`, initial values are used to help the optimizer reaching the optimum parameters. 

I suggest we create a utility function that returns initial values similar to what is done in R. It would only support GEV for now, but could be extended as the need arises. Hopefull scipy will improve its own algorithms in the meantime. 






Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 9b8bd3f8..5abe2ece 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -22,6 +22,7 @@ New features and enhancements
 * New `generic.parametric_quantile` function taking parameters estimated by `generic.fit` as an input.
 * Add support for using probability weighted moments method in `generic.fit` function. Requires the
   `lmoments3` package, which is not included in dependencies because it is unmaintained. Install manually if needed.
+* Implemented `_fit_start` utility function providing initial conditions for statistical distribution parameters estimation, reducing the likelihood of poor fits.
 * Added support for indicators based on hourly (1H) inputs, and a first hourly indicator called `max_pr_intensity`.
   returning hourly precipitation intensity.
 * Indicator instances can be retrieved through their class with the `get_instance()` class method.
diff --git a/xclim/indices/generic.py b/xclim/indices/generic.py
index 6cb65d16..2ee46472 100644
--- a/xclim/indices/generic.py
+++ b/xclim/indices/generic.py
@@ -152,7 +152,8 @@ def fit(da: xr.DataArray, dist: str = "norm", method="ML"):
 
         # Estimate parameters
         if method == "ML":
-            params = dc.fit(x)
+            args, kwargs = _fit_start(x, dist)
+            params = dc.fit(x, *args, **kwargs)
         elif method == "PWM":
             params = list(lm3dc.lmom_fit(x).values())
 
@@ -399,6 +400,42 @@ def get_lm3_dist(dist):
     return getattr(lmoments3.distr, _lm3_dist_map[dist])
 
 
+def _fit_start(x, dist):
+    """Return initial values for distribution parameters.
+
+    Providing the ML fit method initial values can help the optimizer find the global optimum.
+
+    Parameters
+    ----------
+    x : array-like
+      Input data.
+    dist : str
+      Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
+      (see scipy.stats). Only `genextreme` and `weibull_exp` distributions are supported.
+
+    References
+    ----------
+    Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer-Verlag, London, U.K., 208pp
+    Cohen & Whittle, (1988) "Parameter Estimation in Reliability and Life Span Models", p. 25 ff, Marcel Dekker.
+    """
+    x = np.asarray(x)
+    m = x.mean()
+    v = x.var()
+
+    if dist == "genextreme":
+        s = np.sqrt(6 * v) / np.pi
+        return (0.1,), {"loc": m - 0.57722 * s, "scale": s}
+
+    if dist in ("weibull_min"):
+        s = x.std()
+        loc = x.min() - 0.01 * s
+        chat = np.pi / np.sqrt(6) / (np.log(x - loc)).std()
+        scale = ((x - loc) ** chat).mean() ** (1 / chat)
+        return (chat,), {"loc": loc, "scale": scale}
+
+    return (), {}
+
+
 binary_ops = {">": "gt", "<": "lt", ">=": "ge", "<=": "le"}
 
 


Test patch:
diff --git a/tests/test_generic.py b/tests/test_generic.py
index 323b5a45..4a7178be 100644
--- a/tests/test_generic.py
+++ b/tests/test_generic.py
@@ -24,6 +24,80 @@ class TestFA(object):
             coords={"time": time, "x": cx, "y": cy},
         )
 
+        self.weibull_min = xr.DataArray(
+            [
+                4836.6,
+                823.6,
+                3131.7,
+                1343.4,
+                709.7,
+                610.6,
+                3034.2,
+                1973,
+                7358.5,
+                265,
+                4590.5,
+                5440.4,
+                4613.7,
+                4763.1,
+                115.3,
+                5385.1,
+                6398.1,
+                8444.6,
+                2397.1,
+                3259.7,
+                307.5,
+                4607.4,
+                6523.7,
+                600.3,
+                2813.5,
+                6119.8,
+                6438.8,
+                2799.1,
+                2849.8,
+                5309.6,
+                3182.4,
+                705.5,
+                5673.3,
+                2939.9,
+                2631.8,
+                5002.1,
+                1967.3,
+                2810.4,
+                2948,
+                6904.8,
+            ],
+            dims=("time",),
+            coords={"time": xr.IndexVariable("time", np.arange(40))},
+        )
+
+        self.genextreme = xr.DataArray(
+            [
+                279,
+                302,
+                450,
+                272,
+                401,
+                222,
+                311,
+                327,
+                294,
+                299,
+                348,
+                286,
+                492,
+                296,
+                227,
+                437,
+                340,
+                376,
+                444,
+                177,
+            ],
+            dims=("time",),
+            coords={"time": xr.IndexVariable("time", np.arange(20))},
+        )
+
     def test_fit(self):
         p = generic.fit(self.da, "lognorm")
 
@@ -37,6 +111,16 @@ class TestFA(object):
         assert cdf.shape == (self.nx, self.ny)
         assert p.attrs["estimator"] == "Maximum likelihood"
 
+    def test_weibull_min_fit(self):
+        """Check ML fit with a series that leads to poor values without good initial conditions."""
+        p = generic.fit(self.weibull_min, "weibull_min")
+        np.testing.assert_allclose(p, (1.7760067, -322.092552, 4355.262679), 1e-5)
+
+    def test_genextreme_fit(self):
+        """Check ML fit with a series that leads to poor values without good initial conditions."""
+        p = generic.fit(self.genextreme, "genextreme")
+        np.testing.assert_allclose(p, (0.20949, 297.954091, 75.7911863), 1e-5)
+
     def test_fa(self):
         T = 10
         q = generic.fa(self.da, T, "lognorm")
diff --git a/tests/test_units.py b/tests/test_units.py
index 6e71e0a6..bb239a22 100644
--- a/tests/test_units.py
+++ b/tests/test_units.py
@@ -138,3 +138,8 @@ class TestCheckUnits:
 
             with pytest.raises(ValidationError):
                 check_units("m3", "[discharge]")
+
+    @pytest.mark.xfail
+    def test_user_error(self):
+        with pytest.raises(ValidationError):
+            check_units("deg C", "[temperature]")