Fix issue with zero precipitation for distributions that use Gamma function

Author: e-baumer

data/monthly_data_v2.csv

@@ -0,0 +1,194 @@
+"date","precip"
+"1990-01-01",161
+"1990-02-01",40.4
+"1990-03-01",65
+"1990-04-01",83.8
+"1990-05-01",42.8
+"1990-06-01",53.4
+"1990-07-01",43.2
+"1990-08-01",37.4
+"1990-09-01",10.6
+"1990-10-01",1
+"1990-11-01",44.6
+"1990-12-01",94.6
+"1991-01-01",143
+"1991-02-01",25.8
+"1991-03-01",204.4
+"1991-04-01",247.4
+"1991-05-01",33.4
+"1991-06-01",36
+"1991-07-01",20
+"1991-08-01",27.6
+"1991-09-01",8.4
+"1991-10-01",24.6
+"1991-11-01",96.4
+"1991-12-01",192
+"1992-01-01",145.2
+"1992-02-01",38.8
+"1992-03-01",56.6
+"1992-04-01",103.6
+"1992-05-01",49.8
+"1992-06-01",54.2
+"1992-07-01",64.4
+"1992-08-01",0.4
+"1992-09-01",108.2
+"1992-10-01",13.6
+"1992-11-01",44.8
+"1992-12-01",49.2
+"1993-01-01",125
+"1993-02-01",5.2
+"1993-03-01",80.8
+"1993-04-01",113.4
+"1993-05-01",81.6
+"1993-06-01",52
+"1993-07-01",21.4
+"1993-08-01",22.2
+"1993-09-01",19.8
+"1993-10-01",50.8
+"1993-11-01",26.6
+"1993-12-01",141.6
+"1994-01-01",279.2
+"1994-02-01",131.2
+"1994-03-01",103.2
+"1994-04-01",31
+"1994-05-01",18.4
+"1994-06-01",14.4
+"1994-07-01",0.6
+"1994-08-01",17.4
+"1994-09-01",20.8
+"1994-10-01",57.2
+"1994-11-01",161.4
+"1994-12-01",221.8
+"1995-01-01",91
+"1995-02-01",34.8
+"1995-03-01",12.4
+"1995-04-01",48.8
+"1995-05-01",52.8
+"1995-06-01",89.2
+"1995-07-01",9.8
+"1995-08-01",23.4
+"1995-09-01",14.6
+"1995-10-01",129.8
+"1995-11-01",133.2
+"1995-12-01",78.2
+"1996-01-01",174
+"1996-02-01",106.6
+"1996-03-01",78.8
+"1996-04-01",48.4
+"1996-05-01",10
+"1996-06-01",43.8
+"1996-07-01",32.4
+"1996-08-01",55.2
+"1996-09-01",57
+"1996-10-01",21.4
+"1996-11-01",85.4
+"1996-12-01",82.6
+"1997-01-01",99
+"1997-02-01",57.2
+"1997-03-01",20.6
+"1997-04-01",3.2
+"1997-05-01",9
+"1997-06-01",112.6
+"1997-07-01",1
+"1997-08-01",91.2
+"1997-09-01",32.4
+"1997-10-01",55.2
+"1997-11-01",61.8
+"1997-12-01",79.2
+"1998-01-01",184.4
+"1998-02-01",216.8
+"1998-03-01",41.8
+"1998-04-01",16.8
+"1998-05-01",52
+"1998-06-01",121.8
+"1998-07-01",91.4
+"1998-08-01",15.8
+"1998-09-01",55.6
+"1998-10-01",55.2
+"1998-11-01",326.6
+"1998-12-01",62.6
+"1999-01-01",77.4
+"1999-02-01",131.8
+"1999-03-01",47.6
+"1999-04-01",195.2
+"1999-05-01",240.8
+"1999-06-01",87.6
+"1999-07-01",3.6
+"1999-08-01",4
+"1999-09-01",23.4
+"1999-10-01",57.4
+"1999-11-01",31.4
+"1999-12-01",172.2
+"2000-01-01",40.2
+"2000-02-01",272.4
+"2000-03-01",162.2
+"2000-04-01",36.4
+"2000-05-01",61.2
+"2000-06-01",6.2
+"2000-07-01",27.8
+"2000-08-01",104.6
+"2000-09-01",103.6
+"2000-10-01",306.4
+"2000-11-01",57.8
+"2000-12-01",479.8
+"2001-01-01",315.8
+"2001-02-01",157
+"2001-03-01",149.2
+"2001-04-01",100.2
+"2001-05-01",69.2
+"2001-06-01",10.6
+"2001-07-01",12.4
+"2001-08-01",76.8
+"2001-09-01",18.6
+"2001-10-01",117.8
+"2001-11-01",15.6
+"2001-12-01",132.4
+"2002-01-01",343.4
+"2002-02-01",162.2
+"2002-03-01",118.2
+"2002-04-01",149
+"2002-05-01",36.4
+"2002-06-01",122
+"2002-07-01",51.2
+"2002-08-01",0.2
+"2002-09-01",5.6
+"2002-10-01",21.4
+"2002-11-01",116.8
+"2002-12-01",50.8
+"2003-01-01",281.6
+"2003-02-01",250.6
+"2003-03-01",172.4
+"2003-04-01",104.6
+"2003-05-01",44
+"2003-06-01",59.4
+"2003-07-01",39
+"2003-08-01",0.2
+"2003-09-01",13.2
+"2003-10-01",22.6
+"2003-11-01",104.8
+"2003-12-01",19.2
+"2004-01-01",143.8
+"2004-02-01",15.8
+"2004-03-01",164.4
+"2004-04-01",14.2
+"2004-05-01",28.6
+"2004-06-01",16.2
+"2004-07-01",14.2
+"2004-08-01",93.2
+"2004-09-01",27.6
+"2004-10-01",5.8
+"2004-11-01",142
+"2004-12-01",124.4
+"2005-01-01",322.4
+"2005-02-01",274.8
+"2005-03-01",161.4
+"2005-04-01",146.2
+"2005-05-01",248.8
+"2005-06-01",34.6
+"2005-07-01",9.2
+"2005-08-01",27.6
+"2005-09-01",54.2
+"2005-10-01",55.8
+"2005-11-01",74.2
+"2005-12-01",36.2
+"2006-01-01",51.4

docs/wmo_199_en.pdf

@@ -31,6 +31,7 @@ class BaseStandardIndex(object):
 
     def __init__(self):
         self.distrb = None
+        self.non_zero_distr = ['gam', 'pe3']
 
     @staticmethod
     def rolling_window_sum(df: pd.DataFrame, span: int, window_type, center, **kwargs):
@@ -156,11 +157,19 @@ def fit_distribution(self, data: np.array, dist_type: str, fit_type: str='lmom',
         '''
         Fit given distribution to historical precipitation data.
         The fit is accomplished using either L-moments or MLE (Maximum Likelihood Estimation).
+
+        For distributions that use the Gamma Function (Gamma and Pearson 3)
         '''
 
         # Get distribution type
         self.distrb = getattr(distr, dist_type)
 
+        # Determine zeros if distribution can not handle x = 0
+        p_zero = None
+        if dist_type in self.non_zero_distr:
+            p_zero = data[data == 0].shape[0] / data.shape[0]
+            data = data[data != 0]
+
         # Fit distribution
         if fit_type == 'lmom':
             params = self.distrb.lmom_fit(data, **kwargs)
@@ -171,7 +180,7 @@ def fit_distribution(self, data: np.array, dist_type: str, fit_type: str='lmom',
         else:
             raise AttributeError(f"{fit_type} is not an option. Option fit_types are mle and lmom")
 
-        return params
+        return params, p_zero
 
     def calculate(self, df: pd.DataFrame, date_col: str, precip_cols: list, freq: str="M",
                   freq_col: str=None, fit_type: str='lmom', dist_type: str='gam', **dist_kwargs):
@@ -245,6 +254,9 @@ def calculate(self, df: pd.DataFrame, date_col: str, precip_cols: list, freq: st
         if isinstance(precip_cols, str):
             precip_cols = [precip_cols]
 
+        # for p_col in precip_cols:
+            # df.loc[df[p_col] == 0.0, p_col] = 0.001
+
         self._df_copy = df[[date_col] + precip_cols]
         self._df_copy[date_col] = pd.to_datetime(self._df_copy[date_col])
 
@@ -274,12 +286,12 @@ def calculate(self, df: pd.DataFrame, date_col: str, precip_cols: list, freq: st
                 precip_sorted = np.sort(precip_single)[::-1]
 
                 # Fit distribution for particular series and month
-                params = self.fit_distribution(
+                params, p_zero = self.fit_distribution(
                     precip_sorted, dist_type, fit_type, **dist_kwargs
                 )
 
                 # Calculate SPI/SPEI
-                spi = self.cdf_to_ppf(precip_single, params)
+                spi = self.cdf_to_ppf(precip_single, params, p_zero)
                 idx_col = f"{p}_calculated_index"
                 precip_single_df[idx_col] = spi
                 precip_single_df = precip_single_df[[date_col, idx_col]]
@@ -294,7 +306,7 @@ def calculate(self, df: pd.DataFrame, date_col: str, precip_cols: list, freq: st
 
         return df_all
 
-    def cdf_to_ppf(self, data, params):
+    def cdf_to_ppf(self, data, params, p_zero):
         '''
         Take the specific distributions fitted parameters and calculate the
         cdf. Apply the inverse normal distribution to the cdf to get the SPI
@@ -304,7 +316,10 @@ def cdf_to_ppf(self, data, params):
         '''
 
         # Calculate the CDF of observed precipitation on a given time scale
-        cdf = self.distrb.cdf(data, **params)
+        if p_zero:
+            cdf = p_zero + (1 - p_zero) * self.distrb.cdf(data, **params)
+        else:
+            cdf = self.distrb.cdf(data, **params)
 
         # Apply inverse normal distribution
         norm_ppf = scs.norm.ppf(cdf)