Add xclim.indices.stats.fit method: MSE (maximum product of spacings) (…

…#2077)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #2078 
- [ ] Tests for the changes have been added (for bug fixes / features)
- [ ] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGELOG.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

### What kind of change does this PR introduce?

* Adds a fit method: MSE, using
[scipy.stats.fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.fit.html#ra4db2bb0bf1d-1).
Note that this is not available in the `dist` object.

### Does this PR introduce a breaking change?
- No.

### Other information:
- [Maximum spacing
estimation](https://en.wikipedia.org/wiki/Maximum_spacing_estimation#Sensitivity)
- A comparative analysis of L-moments, maximum likelihood, and maximum
product of spacing methods for the four-parameter kappa distribution in
extreme value analysis:
https://www.nature.com/articles/s41598-024-84056-1

CHANGELOG.rst

@@ -30,6 +30,7 @@ New features and enhancements
 * Missing values method "pct" and "at_least_n" now accept a new "subfreq" option that allows to compute the missing mask in two steps. When given, the algorithm is applied at this "subfreq" resampling frequency first and then the result is resampled at the target indicator frequency. In the output, a period is invalid if any of its subgroup where flagged as invalid by the chosen method. (:pull:`2058`, :issue:`1820`).
 * ``scipy.stats.rv_continuous`` instances can now be given directly as the ``dist`` argument in ``standardized_precipitation_index`` and ``standardized_precipitation_evapotranspiration_index`` indicators. This includes `lmoments3` distributions when specifying ``method="PWM"``. (:issue:`2043`, :pull:`2045`).
 * Time selection in ``xclim.core.calendar.select_time`` and the ``**indexer`` argument of indicators now supports day-of-year bounds given as DataArrays with spatial and/or temporal dimensions. (:issue:`1987`, :pull:`2055`).
+* Maximum Spacing Estimation method for distribution fitting has been added to `xclim.indices.stats.fit` (:issue:`2078`, :pull:`2077`)
 
 Bug fixes
 ^^^^^^^^^

src/xclim/indices/stats.py

@@ -55,6 +55,16 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs):
         params = dist.fit(x, *args, method="mle", **kwargs, **fitkwargs)
     elif method == "MM":
         params = dist.fit(x, method="mm", **fitkwargs)
+    elif method in ["MSE", "MPS"]:
+        args, guess = _fit_start(x, dist.name, **fitkwargs)
+        param_info = dist.shapes
+        for i, arg in enumerate(args):
+            guess[param_info[i]] = arg
+
+        fitresult = scipy.stats.fit(
+            dist=dist, data=x, method="mse", guess=guess, **fitkwargs
+        )
+        params = fitresult.params
     elif method == "PWM":
         # lmoments3 will raise an error if only dist.numargs + 2 values are provided
         if len(x) <= dist.numargs + 2:
@@ -107,10 +117,14 @@ def fit(
     dist : str or rv_continuous distribution object
         Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
         (see :py:mod:scipy.stats for full list) or the distribution object itself.
-    method : {"ML", "MLE", "MM", "PWM", "APP"}
-        Fitting method, either maximum likelihood (ML or MLE), method of moments (MM) or approximate method (APP).
+    method : {"ML", "MLE", "MM", "PWM", "APP", "MSE", "MPS"}
+        Fitting method, either maximum likelihood (ML or MLE), method of moments (MM), maximum product of spacings (MSE or MPS)
+        or approximate method (APP).
         Can also be the probability weighted moments (PWM), also called L-Moments, if a compatible `dist` object is passed.
-        The PWM method is usually more robust to outliers.
+        The PWM method is usually more robust to outliers. The MSE method is more consistent than the MLE method, although
+        it can be more sensitive to repeated data.
+        For the MSE method, each variable parameter must be given finite bounds
+        (provided with keyword argument bounds={'param_name':(min,max),...}).
     dim : str
         The dimension upon which to perform the indexing (default: "time").
     **fitkwargs : dict
@@ -131,6 +145,8 @@ def fit(
         "ML": "maximum likelihood",
         "MM": "method of moments",
         "MLE": "maximum likelihood",
+        "MSE": "maximum product of spacings",
+        "MPS": "maximum product of spacings",
         "PWM": "probability weighted moments",
         "APP": "approximative method",
     }

tests/test_stats.py

@@ -2,9 +2,12 @@
 
 from __future__ import annotations
 
+from functools import partial
+
 import numpy as np
 import pytest
 import xarray as xr
+from scipy.optimize import differential_evolution
 from scipy.stats import lognorm, norm
 
 from xclim.indices import stats
@@ -157,6 +160,24 @@ def test_genextreme_fit(genextreme):
     np.testing.assert_allclose(p, (0.20949, 297.954091, 75.7911863), 1e-5)
 
 
+def test_mse_fit(genextreme):
+    """Check MSE fit with a series that leads to poor values without good initial conditions."""
+    # Use fixed-seed rng to remove randomness of differential_evolution
+    # (alternative: change optimizer to non-stochastic gradient-based)
+    # TODO: change when minimum scipy exceeds 1.15.0 to rng=0
+    optimizer = partial(differential_evolution, seed=0)
+    p = stats.fit(
+        genextreme,
+        "genextreme",
+        "MSE",
+        bounds=dict(c=(0, 1), scale=(0, 100), loc=(200, 400)),
+        optimizer=optimizer,
+    )
+    np.testing.assert_allclose(
+        p, (0.18435517630019815, 293.61049928703073, 86.70937297745427), 1e-3
+    )
+
+
 def test_fa(fitda):
     T = 10
     q = stats.fa(fitda, T, "lognorm")