Proper use of loc parameter in gamma,fisk dists (SPI/SPEI)

[coxipi]

Pull Request Checklist:

• This PR addresses an already opened issue (for bug fixes / features)
  • This PR fixes #xyz
• Tests for the changes have been added (for bug fixes / features)
  • (If applicable) Documentation has been added / updated (for bug fixes / features)
• CHANGES.rst has been updated (with summary of main changes)
  • Link to issue (:issue:number) and pull request (:pull:number) has been added

What kind of change does this PR introduce?

• Spei used to rely on an offset to ensure that distribution of values are positive when using distributions such as gamma, fisk. Now, instead, we properly use the loc parameter in those distributions that generalize the distributions and allow negative values. For instance, the gamma disitribution is a 2-parameter function, but in scipy it's generalized to 3 parameters, with a loc parameter. This extends the support of gamma from $x \in [0,\infty]$ to $x \in [\text{loc}, \infty]$.

The code was already compatible with this idea, just had to change _fit_start so it can compute estimates of the parameters in the cases with negatives values.

Does this PR introduce a breaking change?

• Removing the offset from SPEI. Also, the _fit_start can now compute an estimate of loc which can affect what value we obtain when optimizing 3-parameters distributions. Even when the loc is fixed to 0 with floc=0 in fitkwargs, the estimate of fisk was modified. The previous formulation was obtained by applying more or less the idea in the method of moments, but now I obtained an approximation that's more rigorous.

• The case of zero inflated distribution has been changed: Now, it is the calibration data that determines the probability of zero, as it should.

• Method APP can now only be used if fitkwargs['floc'] is given as input by the user.

Other information:

[coxipi]

I started a new branch spi_loc2 because there were many errors hard to track in spi_loc. Opening a new PR, this is the old one: #1653

[coxipi]

I was just trying to tie ends together, and was exploring the difference between fitting with zero inflated distributions vs. normal distributions.

I'm considering the case where the dataset is inflated in some small value that I inject, but I don't do any special treatment for this value. I realize that things work well unless the inflated value is 0:

pr0 = open_dataset("sdba/CanESM2_1950-2100.nc")
# Take the first 30 july 1rst 
pr0  = pr0.isel(location=1).isel(time=slice(179,365*30,365)).pr.values
pr0 = np.where(pr0>0, pr0, 1e-6)
min_prs = [-1e-10,-1e-20, 0, 1e-20, 1e-18, 1e-14, 1e-10]
print(["min_pr", "Pr(X>=min_pr) [fraction]"])
for min_pr in min_prs:
    pr = pr0.copy()
    pr[0:11] = min_pr 
    a, loc,scale = sc.stats.gamma.fit(pr)
    print([min_pr, sc.stats.gamma.cdf(min_pr,a,loc=loc,scale=scale)])

yields:

['min_pr', 'Pr(X>=min_pr) [fraction]']
[-1e-10, 0.07865019433968551]
[-1e-20, 0.07258213369500696]
[0, 5.552055328806407e-20]
[1e-20, 0.07787936999249587]
[1e-18, 0.07488806373878135]
[1e-14, 0.08371030972954428]
[1e-10, 0.08298754087251842]

I inflated the datasets, so we could expect 'Pr(X>=min_pr)' to be about 33%. We have approximate fits, so it's ok to have 7-8%, If I increase the bound in sc.stats.gamma.cdf(min_pr,a,loc=loc,scale=scale) say by adding 1e-15 which is still well below other values that were not injected, then I get indeed ~30%, except for the case of 0 precipitations.

and we can clearly see something is off with the 0-inflated case if we plot the pdfs:

Using a form of jitter_under_thresh, say that I inflate with precipitation values of 0, but then jitter between 10^(-21) and 10^(-20), I get something more in line with other results.

If I inflate less values, 10 or less, it seems the problem goes away (I also added a small value in the cdf computation, but that doesn't explain the change in behaviour as far I could tell):

import numpy as np
import scipy as sc
import matplotlib.pyplot as plt
pr0 = open_dataset("sdba/CanESM2_1950-2100.nc")
# Take the first 30 july 1rst 
pr0  = pr0.isel(location=1).isel(time=slice(179,365*30,365)).pr.values
pr0 = np.where(pr0>0, pr0, 1e-6)
min_prs = [
    0,
    1e-20,
    1e-17,
    1e-16,
    ]
print(["min_pr", "Pr(X>=min_pr)"])
for min_pr in min_prs:
    pr = pr0.copy()
    pr[0:10] = min_pr
    a, loc,scale = sc.stats.gamma.fit(pr)
    print([min_pr, sc.stats.gamma.cdf(min_pr+1e-21,a,loc=loc,scale=scale)])
    x = np.linspace(0, 1e-15, 1000)
    plt.plot(x,sc.stats.gamma.pdf(x, a, loc=loc, scale=scale),label=min_pr)
    plt.yscale("log")
    plt.xlabel("x")
    plt.ylabel("pdf")
plt.legend(title="min_pr")

Going about 11/30 values, then more cases look like the problematic pr==0.

Anyways, it's just worth keeping this in mind. I think the lesson is to be careful when using numbers where the precision is too low, and 0 values are afllicted by this. If we work with the appropriate zero_inflated options or perform some jitter_under_thresh with a very small value, in both cases, this should not be an issue.

[coxipi]

Apparently the PR was never in draft mode? Just wanted to say it's ready for your eyes. ⚠️

[aulemahal]

I like this! thanks a lot !

[aulemahal]

Suggested change

	offset = convert_units_to("1 mm/d", wb, context="hydro")
	offset = convert_units_to("1000 mm/month", wb, context="hydro")

Isn't this the hardcoded offset of climate_indices ?
https://github.com/monocongo/climate_indices/blob/d108eee982abae06a415e888319b8078af868558/src/climate_indices/indices.py#L280C62-L280C69

[coxipi]

yes, climate_indices and xclim do match, but with the current datasets used in our test, we stil have negative values with the climate_indices offset. I wll wait for a next PR where we show many comparisons to many other libraries to sort this out

[coxipi]

but I agree, I didn't have the correct offset here