From b385dba4f103de70277e9a173c63e71a6d7cfae6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 13 Jan 2025 16:13:25 -0500 Subject: [PATCH 01/13] warn fitkwargs & lmoments, shift gamma with floc --- CHANGELOG.rst | 8 ++++++++ src/xclim/indices/stats.py | 15 +++++++++++++++ 2 files changed, 23 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index a4d922613..61da8f36e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -2,6 +2,14 @@ Changelog ========= +v0.55.0 (unpublished) +-------------------- +Contributors to this version: Éric Dupuis (:user:`coxipi`). + +Internal changes +^^^^^^^^^^^^^^^^ +* There is now a warning stating that `fitkwargs` are not employed when using `lmoments3` distribution. One exception is the use of `floc` which is allowed with the gamma distribution. `floc` is used to shift the distribution before computing fitting parameters with the lmoments3 distribution since `loc=0` is always assumed in the library. + v0.54.0 (2024-12-16) -------------------- Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`), Sascha Hofmann (:user:`saschahofmann`). diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 76a018baa..e42e6d770 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -59,6 +59,21 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): # lmoments3 will raise an error if only dist.numargs + 2 values are provided if len(x) <= dist.numargs + 2: return np.asarray([np.nan] * nparams) + if (type(dist).__name__ != "GammaGen" and len(fitkwargs.keys()) != 0) or ( + type(dist).__name__ == "GammaGen" + and set(fitkwargs.keys()) - {"floc"} != set() + ): + warnings.warn( + "Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution." + ) + if "floc" in fitkwargs and type(dist).__name__ == "GammaGen": + # lmoments3 assumes `loc` is 0, so `x` may need to be shifted + # note that `floc` must already be in appropriate units for `x` + params = dist.lmom_fit(x - fitkwargs["floc"]) + params["loc"] = fitkwargs["floc"] + params = list(params.values()) + else: + params = list(dist.lmom_fit(x).values()) params = list(dist.lmom_fit(x).values()) elif method == "APP": args, kwargs = _fit_start(x, dist.name, **fitkwargs) From 6209134b70838215d21d408fe294f18873fe2420 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 13 Jan 2025 16:16:57 -0500 Subject: [PATCH 02/13] update CHANGELOG --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 61da8f36e..b4157122a 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -8,7 +8,7 @@ Contributors to this version: Éric Dupuis (:user:`coxipi`). Internal changes ^^^^^^^^^^^^^^^^ -* There is now a warning stating that `fitkwargs` are not employed when using `lmoments3` distribution. One exception is the use of `floc` which is allowed with the gamma distribution. `floc` is used to shift the distribution before computing fitting parameters with the lmoments3 distribution since `loc=0` is always assumed in the library. +* There is now a warning stating that `fitkwargs` are not employed when using `lmoments3` distribution. One exception is the use of `floc` which is allowed with the gamma distribution. `floc` is used to shift the distribution before computing fitting parameters with the lmoments3 distribution since `loc=0` is always assumed in the library. (:issue:`2043`, :pull:`2045`). v0.54.0 (2024-12-16) -------------------- From 1fb25bd42f00ae9d987824f407733add4db6ea20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Tue, 14 Jan 2025 09:37:34 -0500 Subject: [PATCH 03/13] remove old line Co-authored-by: Pascal Bourgault --- src/xclim/indices/stats.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index e42e6d770..897fa091b 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -74,7 +74,6 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): params = list(params.values()) else: params = list(dist.lmom_fit(x).values()) - params = list(dist.lmom_fit(x).values()) elif method == "APP": args, kwargs = _fit_start(x, dist.name, **fitkwargs) kwargs.setdefault("loc", 0) From 8eb2b51db385d8adc0e6aaa639586756ea9a4fd3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 09:46:43 -0500 Subject: [PATCH 04/13] fitkwargs&PWM=error, allow PWM in SPI/SPEI --- src/xclim/indices/_agro.py | 6 ++++-- src/xclim/indices/stats.py | 5 +++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 52726c3f9..baa08662a 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1150,6 +1150,7 @@ def standardized_precipitation_index( uses a deterministic function that does not involve any optimization. fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). + If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed. cal_start : DateStr, optional Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. Default option `None` means that the calibration period begins at the start of the input dataset. @@ -1281,11 +1282,12 @@ def standardized_precipitation_evapotranspiration_index( resampling, the window is an integer number of months. dist : {'gamma', 'fisk'} Name of the univariate distribution. (see :py:mod:`scipy.stats`). - method : {'APP', 'ML'} + method : {'APP', 'ML', 'PWM'} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that doesn't involve any optimization. fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). + If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed. cal_start : DateStr, optional Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. Default option `None` means that the calibration period begins at the start of the input dataset. @@ -1311,7 +1313,7 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP", "PWM"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 897fa091b..b021857e7 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -63,7 +63,7 @@ def _fitfunc_1d(arr, *, dist, nparams, method, **fitkwargs): type(dist).__name__ == "GammaGen" and set(fitkwargs.keys()) - {"floc"} != set() ): - warnings.warn( + raise ValueError( "Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution." ) if "floc" in fitkwargs and type(dist).__name__ == "GammaGen": @@ -884,8 +884,9 @@ def standardized_index( The approximate method uses a deterministic function that doesn't involve any optimization. zero_inflated : bool If True, the zeroes of `da` are treated separately. - fitkwargs : dict + fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). + If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed. cal_start : DateStr, optional Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. Default option `None` means that the calibration period begins at the start of the input dataset. From 2ef071c03d07ffc4670d427326a816c6340e4638 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 15:58:21 -0500 Subject: [PATCH 05/13] add tests --- src/xclim/indices/_agro.py | 6 +-- src/xclim/indices/stats.py | 2 +- tests/test_indices.py | 75 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 79 insertions(+), 4 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index baa08662a..fd9355959 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1145,7 +1145,7 @@ def standardized_precipitation_index( i.e. a monthly resampling, the window is an integer number of months. dist : {"gamma", "fisk"} Name of the univariate distribution. (see :py:mod:`scipy.stats`). - method : {"APP", "ML"} + method : {"APP", "ML", "PWM"} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that does not involve any optimization. fitkwargs : dict, optional @@ -1219,7 +1219,7 @@ def standardized_precipitation_index( >>> spi_3_fitted = standardized_precipitation_index(pr, params=params) """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( @@ -1313,7 +1313,7 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP", "PWM"]} + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index b021857e7..477f8d7cf 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -808,7 +808,7 @@ def standardized_index_fit_params( "Pass a value for `floc` in `fitkwargs`." ) - dist_and_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} + dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} dist = get_dist(dist) if dist.name not in dist_and_methods: raise NotImplementedError(f"The distribution `{dist.name}` is not supported.") diff --git a/tests/test_indices.py b/tests/test_indices.py index 22978b49d..c8444fec4 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -660,6 +660,38 @@ class TestStandardizedIndices: [-1.08627775, -0.46491398, -0.77806462, 0.31759127, 0.03794528], 2e-2, ), + ( + "D", + 1, + "gamma", + "PWM", + [-0.13002, 1.346689, 0.965731, 0.245408, -0.427896], + 2e-2, + ), + ( + "D", + 12, + "gamma", + "PWM", + [-0.209411, -0.086357, 0.636851, 1.022608, 0.634409], + 2e-2, + ), + ( + "MS", + 1, + "gamma", + "PWM", + [1.364243, 1.478565, 1.915559, -3.055828, 0.905304], + 2e-2, + ), + ( + "MS", + 12, + "gamma", + "PWM", + [0.577214, 1.522867, 1.634222, 0.967847, 0.689001], + 2e-2, + ), ], ) def test_standardized_precipitation_index( @@ -671,6 +703,13 @@ def test_standardized_precipitation_index( and Version(__numpy_version__) < Version("2.0.0") ): pytest.skip("Skipping SPI/ML/D on older numpy") + + # change `dist` to a lmoments3 object if needed + if method == "PWM": + lmom = pytest.importorskip("lmoments3.distr") + scipy2lmom = {"gamma": "gam"} + dist = getattr(lmom, scipy2lmom[dist]) + ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1) if freq == "D": # to compare with ``climate_indices`` @@ -922,6 +961,42 @@ def test_zero_inflated(self, open_dataset): np.all(np.not_equal(spid[False].values, spid[True].values)), True ) + def test_PWM_and_fitkwargs(self, open_dataset): + pr = ( + open_dataset("sdba/CanESM2_1950-2100.nc") + .isel(location=1) + .sel(time=slice("1950", "1980")) + ).pr + + lmom = pytest.importorskip("lmoments3.distr") + # for now, only one function used + scipy2lmom = {"gamma": "gam"} + dist = getattr(lmom, scipy2lmom["gamma"]) + fitkwargs = {"floc": 0} + input_params = dict( + freq=None, + window=1, + method="PWM", + dist=dist, + fitkwargs=fitkwargs, + # doy_bounds=(180, 180), + ) + # this should not cause a problem + params_d0 = xci.stats.standardized_index_fit_params(pr, **input_params).isel( + dayofyear=0 + ) + np.testing.assert_allclose( + params_d0, np.array([5.63e-01, 0, 3.37e-05]), rtol=0, atol=2e-2 + ) + # this should cause a problem + fitkwargs["fscale"] = 1 + input_params["fitkwargs"] = fitkwargs + with pytest.raises( + ValueError, + match="Lmoments3 does not use `fitkwargs` arguments, except for `floc` with the Gamma distribution.", + ): + xci.stats.standardized_index_fit_params(pr, **input_params) + class TestDailyFreezeThawCycles: @pytest.mark.parametrize( From dbce9669227f78a032d53e09781197d9c7ebb091 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 14 Jan 2025 16:05:20 -0500 Subject: [PATCH 06/13] update CHANGELOG --- CHANGELOG.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index e01f7af36..ea7a24d95 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -13,6 +13,7 @@ New features and enhancements Internal changes ^^^^^^^^^^^^^^^^ * There is now a warning stating that `fitkwargs` are not employed when using `lmoments3` distribution. One exception is the use of `floc` which is allowed with the gamma distribution. `floc` is used to shift the distribution before computing fitting parameters with the lmoments3 distribution since `loc=0` is always assumed in the library. (:issue:`2043`, :pull:`2045`). +* `PWM` method (`lmoments3`) is now available to be used with the `gamma` distribution in ``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`2043`, :pull:`2045`). v0.54.0 (2024-12-16) -------------------- From c5d02593ec46de3633d9b02cab97dd4c0dbdc955 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Wed, 15 Jan 2025 13:18:00 -0500 Subject: [PATCH 07/13] remove conditions on PWM --- src/xclim/indices/_agro.py | 14 ++++++++------ src/xclim/indices/stats.py | 17 ++++++++++------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index fd9355959..9637a1d50 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1143,8 +1143,8 @@ def standardized_precipitation_index( window : int Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, i.e. a monthly resampling, the window is an integer number of months. - dist : {"gamma", "fisk"} - Name of the univariate distribution. (see :py:mod:`scipy.stats`). + dist : {'gamma', 'fisk'} + Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. method : {"APP", "ML", "PWM"} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that does not involve any optimization. @@ -1225,7 +1225,8 @@ def standardized_precipitation_index( raise NotImplementedError( f"{method} method is not implemented for {dist} distribution." ) - else: + # Constraints on distributions except for PWM + elif method != "PWM": raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Precipitation is expected to be zero-inflated @@ -1281,7 +1282,7 @@ def standardized_precipitation_evapotranspiration_index( Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, i.e. a monthly resampling, the window is an integer number of months. dist : {'gamma', 'fisk'} - Name of the univariate distribution. (see :py:mod:`scipy.stats`). + Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. method : {'APP', 'ML', 'PWM'} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that doesn't involve any optimization. @@ -1313,13 +1314,14 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( f"{method} method is not implemented for {dist} distribution" ) - else: + # Constraints on distributions except for PWM + elif method != "PWM": raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Water budget is not expected to be zero-inflated diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 477f8d7cf..e2c222ad5 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -808,14 +808,17 @@ def standardized_index_fit_params( "Pass a value for `floc` in `fitkwargs`." ) - dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_and_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} dist = get_dist(dist) - if dist.name not in dist_and_methods: - raise NotImplementedError(f"The distribution `{dist.name}` is not supported.") - if method not in dist_and_methods[dist.name]: - raise NotImplementedError( - f"The method `{method}` is not supported for distribution `{dist.name}`." - ) + if method != "PWM": + if dist.name not in dist_and_methods: + raise NotImplementedError( + f"The distribution `{dist.name}` is not supported." + ) + if method not in dist_and_methods[dist.name]: + raise NotImplementedError( + f"The method `{method}` is not supported for distribution `{dist.name}`." + ) da, group = preprocess_standardized_index(da, freq, window, **indexer) if zero_inflated: prob_of_zero = da.groupby(group).map( From 6052fe553c25254144545eef14a127d889cc55b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 20 Jan 2025 10:57:13 -0500 Subject: [PATCH 08/13] no PWM constraint for SPI either --- src/xclim/indices/_agro.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 9637a1d50..0601c3f0c 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1219,11 +1219,12 @@ def standardized_precipitation_index( >>> spi_3_fitted = standardized_precipitation_index(pr, params=params) """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + + dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( - f"{method} method is not implemented for {dist} distribution." + f"{method} method is not implemented for {dist} distribution" ) # Constraints on distributions except for PWM elif method != "PWM": From 7133fcc1c4a42a626e75a421e48da6f648f4a6e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:14:18 -0500 Subject: [PATCH 09/13] no constraints on any rv_continuous function --- src/xclim/indices/_agro.py | 51 +++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index 0601c3f0c..c08d138e9 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -6,6 +6,7 @@ import numpy as np import xarray +from scipy.stats import rv_continuous import xclim.indices.run_length as rl from xclim.core import DateStr, DayOfYearStr, Quantified @@ -1122,7 +1123,7 @@ def standardized_precipitation_index( pr: xarray.DataArray, freq: str | None = "MS", window: int = 1, - dist: str = "gamma", + dist: str | rv_continuous = "gamma", method: str = "ML", fitkwargs: dict | None = None, cal_start: DateStr | None = None, @@ -1143,11 +1144,11 @@ def standardized_precipitation_index( window : int Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, i.e. a monthly resampling, the window is an integer number of months. - dist : {'gamma', 'fisk'} - Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. + dist : {'gamma', 'fisk'} or `rv_continuous` function + Name of the univariate distribution, or a callable `rv_continuous` (see :py:mod:`scipy.stats`). method : {"APP", "ML", "PWM"} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method - uses a deterministic function that does not involve any optimization. + uses a deterministic function that does not involve any optimization. `PWM` should be used with a `lmoments3` distribution. fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed. @@ -1221,14 +1222,14 @@ def standardized_precipitation_index( fitkwargs = fitkwargs or {} dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} - if dist in dist_methods: - if method not in dist_methods[dist]: - raise NotImplementedError( - f"{method} method is not implemented for {dist} distribution" - ) - # Constraints on distributions except for PWM - elif method != "PWM": - raise NotImplementedError(f"{dist} distribution is not yet implemented.") + if isinstance(dist, str): + if dist in dist_methods: + if method not in dist_methods[dist]: + raise NotImplementedError( + f"{method} method is not implemented for {dist} distribution" + ) + else: + raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Precipitation is expected to be zero-inflated zero_inflated = True @@ -1257,7 +1258,7 @@ def standardized_precipitation_evapotranspiration_index( wb: xarray.DataArray, freq: str | None = "MS", window: int = 1, - dist: str = "gamma", + dist: str | rv_continuous = "gamma", method: str = "ML", fitkwargs: dict | None = None, cal_start: DateStr | None = None, @@ -1282,11 +1283,11 @@ def standardized_precipitation_evapotranspiration_index( window : int Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, i.e. a monthly resampling, the window is an integer number of months. - dist : {'gamma', 'fisk'} - Name of the univariate distribution. (see :py:mod:`scipy.stats`). All possible distributions are allowed with 'PWM'. - method : {'APP', 'ML', 'PWM'} + dist : {'gamma', 'fisk'} or `rv_continuous` function + Name of the univariate distribution, or a callable `rv_continuous` (see :py:mod:`scipy.stats`). + method : {"APP", "ML", "PWM"} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method - uses a deterministic function that doesn't involve any optimization. + uses a deterministic function that does not involve any optimization. `PWM` should be used with a `lmoments3` distribution. fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). If method is `PWM`, `fitkwargs` should be empty, except for `floc` with `dist`=`gamma` which is allowed. @@ -1316,14 +1317,14 @@ def standardized_precipitation_evapotranspiration_index( fitkwargs = fitkwargs or {} dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} - if dist in dist_methods: - if method not in dist_methods[dist]: - raise NotImplementedError( - f"{method} method is not implemented for {dist} distribution" - ) - # Constraints on distributions except for PWM - elif method != "PWM": - raise NotImplementedError(f"{dist} distribution is not yet implemented.") + if isinstance(dist, str): + if dist in dist_methods: + if method not in dist_methods[dist]: + raise NotImplementedError( + f"{method} method is not implemented for {dist} distribution" + ) + else: + raise NotImplementedError(f"{dist} distribution is not yet implemented.") # Water budget is not expected to be zero-inflated zero_inflated = False From 3f0f24b8f3a93c4e8164590ebb2c01ffcc3d76b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:37:12 -0500 Subject: [PATCH 10/13] update CHANGELOG --- CHANGELOG.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 6a1df537e..7c66248ce 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,6 +23,8 @@ New features and enhancements * `xclim` now tracks energy usage and carbon emissions ("last run", "average", and "total") during CI workflows using the `eco-ci-energy-estimation` GitHub Action. (:pull:`2046`). * ``xclim.testing.helpers.test_timeseries`` now accepts a `calendar` argument that is forwarded to ``xr.cftime_range``. (:pull:`2019`). * New ``xclim.indices.fao_allen98``, exporting the FAO-56 Penman-Monteith equation for potential evapotranspiration (:issue:`2004`, :pull:`2067`). +* There is now a warning stating that `fitkwargs` are not employed when using the `lmoments3` distribution. One exception is the use of `'floc'` which is allowed with the gamma distribution. `'floc'` is used to shift the distribution before computing fitting parameters with the `lmoments3` distribution since ``loc=0`` is always assumed in the library. (:issue:`2043`, :pull:`2045`). +* `rv_continuous` functions can now be given directly as the `dist` argument in ``standardized_precipitation_index`` and ``standardized_precipitation_evapotranspiration_index`` indicators. This includes `lmoments3` distribution when specifying `method="PWM"`. (:issue:`2043`, :pull:`2045`). Internal changes ^^^^^^^^^^^^^^^^ @@ -30,7 +32,6 @@ Internal changes * Adjusted the ``TestOfficialYaml`` test to use a dynamic method for finding the installed location of `xclim`. (:pull:`2028`). * Adjusted two tests for better handling when running in Windows environments. (:pull:`2057`). * There is now a warning stating that `fitkwargs` are not employed when using the `lmoments3` distribution. One exception is the use of `'floc'` which is allowed with the gamma distribution. `'floc'` is used to shift the distribution before computing fitting parameters with the `lmoments3` distribution since ``loc=0`` is always assumed in the library. (:issue:`2043`, :pull:`2045`). -* The `PWM` method (from `lmoments3`) is now available to be used with the `gamma` distribution in ``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`2043`, :pull:`2045`). Bug fixes ^^^^^^^^^ From daf85192df42cfc600a2245d1d9a7ef790f50cf5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Tue, 11 Feb 2025 16:41:19 -0500 Subject: [PATCH 11/13] test rv_continuous --- tests/test_indices.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/test_indices.py b/tests/test_indices.py index b6c0f692e..6de6fc626 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -22,6 +22,7 @@ from numpy import __version__ as __numpy_version__ from packaging.version import Version from pint import __version__ as __pint_version__ +from scipy import stats from xclim import indices as xci from xclim.core import ValidationError @@ -749,6 +750,34 @@ def test_standardized_precipitation_index( np.testing.assert_allclose(spi.values, values, rtol=0, atol=diff_tol) + @pytest.mark.slow + @pytest.mark.parametrize("dist", ["gamma", "fisk"]) + def test_str_vs_rv_continuous(self, open_dataset, dist): + ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1) + window = 1 + method = "ML" + freq = "MS" + + pr = ds.pr.sel(time=slice("1998", "2000")) + pr_cal = ds.pr.sel(time=slice("1950", "1980")) + fitkwargs = {} + + out = [] + for dist0 in [dist, getattr(stats, dist)]: + params = xci.stats.standardized_index_fit_params( + pr_cal, + freq=freq, + window=window, + dist=dist0, + method=method, + fitkwargs=fitkwargs, + zero_inflated=True, + ) + spi = xci.standardized_precipitation_index(pr, params=params) + # Only a few moments before year 2000 are tested + out.append(spi.isel(time=slice(-11, -1, 2))) + assert all(out[0] == out[1]) + # See SPI version @pytest.mark.slow @pytest.mark.parametrize( From 6a8a5a12a6b0d3821ec04a1a716026246d9180fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Wed, 12 Feb 2025 13:48:59 -0500 Subject: [PATCH 12/13] better doc --- src/xclim/indices/_agro.py | 16 +++++++++++----- src/xclim/indices/stats.py | 17 +++++++++++++++-- 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/src/xclim/indices/_agro.py b/src/xclim/indices/_agro.py index cfa49faa1..e155f0c9b 100644 --- a/src/xclim/indices/_agro.py +++ b/src/xclim/indices/_agro.py @@ -1171,18 +1171,22 @@ def standardized_precipitation_index( xarray.DataArray, [unitless] Standardized Precipitation Index. + See Also + -------- + xclim.indices.stats.standardized_index : Standardized Index. + xclim.indices.stats.standardized_index_fit_params : Standardized Index Fit Params. + Notes ----- * N-month SPI / N-day SPI is determined by choosing the `window = N` and the appropriate frequency `freq`. * Supported statistical distributions are: ["gamma", "fisk"], where "fisk" is scipy's implementation of - a log-logistic distribution + a log-logistic distribution * Supported frequencies are daily ("D"), weekly ("W"), and monthly ("MS"). - Weekly frequency will only work if the input array has a "standard" (non-cftime) calendar. + * Weekly frequency will only work if the input array has a "standard" (non-cftime) calendar. * If `params` is given as input, it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. * "APP" method only supports two-parameter distributions. Parameter `loc` needs to be fixed to use method `APP`. - * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in - the inversion to the normal distribution. - * The results from `climate_indices` library can be reproduced with `method = "APP"` and `fitwkargs = {"floc": 0}` + * The results from `climate_indices` library can be reproduced with `method = "APP"` and `fitwkargs = {"floc": 0}`, except for the maximum + and minimum values allowed which are greater in xclim ±8.21, . See `xclim.indices.stats.standardized_index` References ---------- @@ -1313,6 +1317,8 @@ def standardized_precipitation_evapotranspiration_index( See Also -------- standardized_precipitation_index : Standardized Precipitation Index. + xclim.indices.stats.standardized_index : Standardized Index. + xclim.indices.stats.standardized_index_fit_params : Standardized Index Fit Params. """ fitkwargs = fitkwargs or {} diff --git a/src/xclim/indices/stats.py b/src/xclim/indices/stats.py index 085ff7e46..ca2b2874a 100644 --- a/src/xclim/indices/stats.py +++ b/src/xclim/indices/stats.py @@ -788,9 +788,11 @@ def standardized_index_fit_params( Notes ----- Supported combinations of `dist` and `method` are: - * Gamma ("gamma") : "ML", "APP", "PWM" + * Gamma ("gamma") : "ML", "APP" * Log-logistic ("fisk") : "ML", "APP" * "APP" method only supports two-parameter distributions. Parameter `loc` will be set to 0 (setting `floc=0` in `fitkwargs`). + * Otherwise, generic `rv_continuous` methods can be used. This includes distributions from `lmoments3` which should be used with + `method="PWM"`. When using the zero inflated option, : A probability density function :math:`\texttt{pdf}_0(X)` is fitted for :math:`X \neq 0` and a supplementary parameter :math:`\pi` takes into account the probability of :math:`X = 0`. The full probability density @@ -909,11 +911,22 @@ def standardized_index( xarray.DataArray, [unitless] Standardized Precipitation Index. + See Also + -------- + standardized_index_fit_params : Standardized Index Fit Params. + Notes ----- * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in the inversion to the normal distribution. - * ``window``, ``dist``, ``method``, ``zero_inflated`` are only optional if ``params`` is given. + * ``window``, ``dist``, ``method``, ``zero_inflated`` are only optional if ``params`` is given. If `params` is given as input, + it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. + * Supported combinations of `dist` and `method` are: + * Gamma ("gamma") : "ML", "APP" + * Log-logistic ("fisk") : "ML", "APP" + * "APP" method only supports two-parameter distributions. Parameter `loc` will be set to 0 (setting `floc=0` in `fitkwargs`). + * Otherwise, generic `rv_continuous` methods can be used. This includes distributions from `lmoments3` which should be used with + `method="PWM"`. References ---------- From 9c608fff25ecd1261aed4addfc5078eb86f9d97b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 13 Feb 2025 00:25:10 -0500 Subject: [PATCH 13/13] remove comment --- tests/test_indices.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 6de6fc626..835582407 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -1009,7 +1009,6 @@ def test_PWM_and_fitkwargs(self, open_dataset): method="PWM", dist=dist, fitkwargs=fitkwargs, - # doy_bounds=(180, 180), ) # this should not cause a problem params_d0 = xci.stats.standardized_index_fit_params(pr, **input_params).isel(