Added loc scale

lmoments3/distr.py

@@ -26,6 +26,14 @@
 import lmoments3 as lm
 
 
+try:
+    # Scipy >= 1.9
+    from scipy.stats._distn_infrastructure import rv_continuous_frozen
+except ImportError:
+    # Scipy < 1.9
+    from scipy.stats._distn_infrastructure import rv_frozen as rv_continuous_frozen
+
+
 class LmomDistrMixin(object):
     """
     Mixin class to add L-moment methods to :class:`scipy.stats.rv_continous` distribution functions. Distributions using
@@ -123,7 +131,7 @@ def freeze(self, *args, **kwds):
         return LmomFrozenDistr(self, *args, **kwds)
 
 
-class LmomFrozenDistr(scipy.stats.distributions.rv_frozen):
+class LmomFrozenDistr(rv_continuous_frozen):
     """
     Frozen version of the distribution returned by :class:`LmomDistrMixin`. Simply provides additional methods supported
     by the mixin.
@@ -460,7 +468,7 @@ def _lmom_fit(self, lmom_ratios):
         T3 = lmom_ratios[2]
         T4 = lmom_ratios[3]
         if lmom_ratios[1] <= 0 or abs(T3) >= 1 or abs(T4) >= 1 or T4 <= (5 * T3 * T3 - 1) / 4 or \
-                        T4 >= (5 * T3 * T3 + 1) / 6:
+                T4 >= (5 * T3 * T3 + 1) / 6:
             raise ValueError("L-Moments invalid")
 
         G = (1 - 3 * T3) / (1 + T3)
@@ -683,35 +691,37 @@ class WakebyGen(LmomDistrMixin, scipy.stats.rv_continuous):
 
     """
 
-    def _argcheck(self, b, c, d):
-        b = np.asarray(b)
-        c = np.asarray(c)
-        d = np.asarray(d)
-        check = np.where(b + d > 0,
-                         np.where(c == 0, d == 0, True),
-                         (b == c) & (c == d) & (d == 0))
-        np.putmask(check, c > 0, d > 0)
-        np.putmask(check, c < 0, False)
+    def _argcheck(self, b, c, d, a):
+        alpha = np.asarray(a)
+        beta = np.asarray(b)
+        gamma = np.asarray(c)
+        delta = np.asarray(d)
+        cond1 = alpha + beta >= 0
+        cond2 = (alpha + beta) > 0 or (b == 0 and gamma == 0 and delta == 0)
+        cond3 = np.where(gamma == 0, delta == 0, True)
+        cond4 = gamma >= 0
+        cond5 = alpha + gamma >= 0
+        check = cond1 and cond2 and cond3 and cond4 and cond5
         return check
 
-    def _ppf(self, q, b, c, d):
+    def _ppf(self, q, b, c, d, a):
         z = -np.log(1. - q)
         u = np.where(b == 0, z, (1. - np.exp(-b * z)) / b)
         v = np.where(d == 0, z, (1. - np.exp(d * z)) / d)
-        return u - c * v
+        return a * u - c * v
 
-    def _cdf(self, x, b, c, d):
+    def _cdf(self, x, b, c, d, a):
         if hasattr(x, '__iter__'):
-            if hasattr(b, '__iter__'):
+            if hasattr(b, '__iter__') and (np.size(b) == np.size(x)):
                 # Assume x, b, c, d are arrays with matching length
                 result = np.array([self._cdfwak(_, parameters)
-                                   for (_, parameters) in zip(x, zip(b, c, d))])
+                                   for (_, parameters) in zip(x, zip(b, c, d, a))])
             else:
                 # Only x is an array, paras are scalars
-                result = np.array([self._cdfwak(_, [b, c, d])
+                result = np.array([self._cdfwak(_, [b, c, d, a])
                                    for _ in x])
         else:
-            result = self._cdfwak(x, (b, c, d))
+            result = self._cdfwak(x, (b, c, d, a))
         return result
 
     def _cdfwak(self, x, para):
@@ -722,9 +732,8 @@ def _cdfwak(self, x, para):
         ZINCMX = 3
         ZMULT = 0.2
         UFL = -170
-        XI = 0  # stats.rv_continuous deals with scaling
-        A = 1  # stats.rv_continuous deals with scaling
-        B, C, D = para
+        XI = 0  # stats.rv_continuous deals with location
+        B, C, D, A = para
 
         CDFWAK = 0
         if x <= XI:
@@ -810,9 +819,9 @@ def _cdfwak(self, x, para):
                         CDFWAK = 1 - math.exp(-Z)
                     return CDFWAK
 
-    def _pdf(self, x, b, c, d):
-        t = (1. - self._cdf(x, b, c, d))
-        f = t ** (d + 1) / (t ** (b + d) + c)
+    def _pdf(self, x, b, c, d, a):
+        t = (1. - self._cdf(x, b, c, d, a))
+        f = t ** (d + 1) / (a * t ** (b + d) + c)
         return f
 
     def _lmom_fit(self, lmom_ratios):
@@ -864,21 +873,23 @@ def _lmom_fit(self, lmom_ratios):
         para = OrderedDict([('beta', B),
                             ('gamma', C),
                             ('delta', D),
+                            ('alpha', A),
                             ('loc', XI),
-                            ('scale', A)])
+                            ('scale', 1)]
+                           )
         return para
 
-    def _lmom_ratios(self, beta, gamma, delta, loc, scale, nmom):
+    def _lmom_ratios(self, beta, gamma, delta, alpha, loc, scale, nmom):
         if (delta >= 1) \
                 or (beta + delta <= 0 and (beta != 0 or gamma != 0 or delta != 0)) \
-                or (scale == 0 and beta != 0) \
+                or (alpha == 0 and beta != 0) \
                 or (gamma == 0 and delta != 0) \
                 or (gamma < 0) \
-                or (scale + gamma < 0) \
-                or (scale == 0 and gamma == 0):
+                or (alpha + gamma < 0) \
+                or (alpha == 0 and gamma == 0):
             raise ValueError("Invalid parameters")
 
-        Y = scale / (1 + beta)
+        Y = alpha / (1 + beta)
         Z = gamma / (1 - delta)
         xmom = [loc + Y + Z]
         if nmom == 1:
@@ -898,7 +909,7 @@ def _lmom_ratios(self, beta, gamma, delta, loc, scale, nmom):
         return xmom
 
 
-wak = WakebyGen(name='wakeby', shapes='beta, gamma, delta')
+wak = WakebyGen(name='wakeby', shapes='beta, gamma, delta, alpha')
 
 """
 The following distributions are available in `scipy.stats` and are redefined here with an `LmomDistrMixin` to extend the

lmoments3/stats.py

@@ -31,7 +31,7 @@ def AIC(data, distr_name, distr_paras):
 
 
 def AICc(data, distr_name, distr_paras):
-    distr_f = getattr(distr, distr_name.lower())  # scipy rv_continous class
+    distr_f = getattr(distr, distr_name.lower())  # scipy rv_continuous class
 
     AICbase = AIC(data, distr_name, distr_paras)
     k = distr_f.numargs + 2  # Include location and scale in addition to shape parameters
@@ -41,7 +41,7 @@ def AICc(data, distr_name, distr_paras):
 
 
 def BIC(data, distr_name, distr_paras):
-    distr_f = getattr(distr, distr_name.lower())  # scipy rv_continous class
+    distr_f = getattr(distr, distr_name.lower())  # scipy rv_continuous class
 
     NLL = distr_f.nnlf(data, **distr_paras)
     k = distr_f.numargs + 2  # Include location and scale in addition to shape parameters

lmoments3/tests/test_lmoments.py

@@ -150,7 +150,7 @@ class TestPe3(DistributionTestCase):
 
 class TestWak(DistributionTestCase):
     dist = 'wak'
-    paras = {'loc': 0.7928727, 'scale': 2.7855796, 'beta': 0.14, 'gamma': 0, 'delta': 0}
+    paras = {'loc': 0.7928727, 'scale': 1, 'alpha': 2.7855796, 'beta': 0.14, 'gamma': 0, 'delta': 0}
     correct_fit = [0.7928727, 2.7855796, 0.1400000, 0.0000000, 0.0000000]
     correct_qua = [1.404848, 2.632964, 4.806899]
     correct_lmr = [3.2363636, 1.1418182, 0.2738854, 0.1230499]