ENH: stats: Use explicit formulas for sf, logsf, etc in weibull distr…

…ibutions. (scipy#6362)

The functions _sf and _logsf are added to the frechet_r/weibull_min
distribution.

The functions _logpdf, _logcdf and _sf are added to the frechet_l/weibull_max
distribution.

The main advantage of these additions is these distributions will now
give more accurate values for the survival function where the default
formula of 1 - CDF gives 0.

scipy/stats/_continuous_distns.py

@@ -1557,6 +1557,12 @@ def _logpdf(self, x, c):
     def _cdf(self, x, c):
         return -special.expm1(-pow(x, c))
 
+    def _sf(self, x, c):
+        return exp(-pow(x, c))
+
+    def _logsf(self, x, c):
+        return -pow(x, c)
+
     def _ppf(self, q, c):
         return pow(-special.log1p(-q), 1.0/c)
 
@@ -1565,6 +1571,7 @@ def _munp(self, n, c):
 
     def _entropy(self, c):
         return -_EULER / c - log(c) + _EULER + 1
+
 frechet_r = frechet_r_gen(a=0.0, name='frechet_r')
 weibull_min = frechet_r_gen(a=0.0, name='weibull_min')
 
@@ -1597,9 +1604,18 @@ class frechet_l_gen(rv_continuous):
     def _pdf(self, x, c):
         return c*pow(-x, c-1)*exp(-pow(-x, c))
 
+    def _logpdf(self, x, c):
+        return log(c) + special.xlogy(c-1, -x) - pow(-x, c)
+
     def _cdf(self, x, c):
         return exp(-pow(-x, c))
 
+    def _logcdf(self, x, c):
+        return -pow(-x, c)
+
+    def _sf(self, x, c):
+        return -special.expm1(-pow(-x, c))
+
     def _ppf(self, q, c):
         return -pow(-log(q), 1.0/c)
 
@@ -1613,6 +1629,7 @@ def _munp(self, n, c):
 
     def _entropy(self, c):
         return -_EULER / c - log(c) + _EULER + 1
+
 frechet_l = frechet_l_gen(b=0.0, name='frechet_l')
 weibull_max = frechet_l_gen(b=0.0, name='weibull_max')
 

scipy/stats/tests/test_distributions.py

@@ -2082,12 +2082,97 @@ def test_a_is_1_c_is_1(self):
         assert_allclose(logp, expected)
 
 
-class TestWeibullMin(TestCase):
-    # gh-6217
+class TestWeibull(TestCase):
+
     def test_logpdf(self):
+        # gh-6217
         y = stats.weibull_min.logpdf(0, 1)
         assert_equal(y, 0)
 
+    def test_with_maxima_distrib(self):
+        # Tests for weibull_min and weibull_max.
+        # The expected values were computed using the symbolic algebra
+        # program 'maxima' with the package 'distrib', which has
+        # 'pdf_weibull' and 'cdf_weibull'.  The mapping between the
+        # scipy and maxima functions is as follows:
+        # -----------------------------------------------------------------
+        # scipy                              maxima
+        # ---------------------------------  ------------------------------
+        # weibull_min.pdf(x, a, scale=b)     pdf_weibull(x, a, b)
+        # weibull_min.logpdf(x, a, scale=b)  log(pdf_weibull(x, a, b))
+        # weibull_min.cdf(x, a, scale=b)     cdf_weibull(x, a, b)
+        # weibull_min.logcdf(x, a, scale=b)  log(cdf_weibull(x, a, b))
+        # weibull_min.sf(x, a, scale=b)      1 - cdf_weibull(x, a, b)
+        # weibull_min.logsf(x, a, scale=b)   log(1 - cdf_weibull(x, a, b))
+        #
+        # weibull_max.pdf(x, a, scale=b)     pdf_weibull(-x, a, b)
+        # weibull_max.logpdf(x, a, scale=b)  log(pdf_weibull(-x, a, b))
+        # weibull_max.cdf(x, a, scale=b)     1 - cdf_weibull(-x, a, b)
+        # weibull_max.logcdf(x, a, scale=b)  log(1 - cdf_weibull(-x, a, b))
+        # weibull_max.sf(x, a, scale=b)      cdf_weibull(-x, a, b)
+        # weibull_max.logsf(x, a, scale=b)   log(cdf_weibull(-x, a, b))
+        # -----------------------------------------------------------------
+        x = 1.5
+        a = 2.0
+        b = 3.0
+
+        # weibull_min
+
+        p = stats.weibull_min.pdf(x, a, scale=b)
+        assert_allclose(p, np.exp(-0.25)/3)
+
+        lp = stats.weibull_min.logpdf(x, a, scale=b)
+        assert_allclose(lp, -0.25 - np.log(3))
+
+        c = stats.weibull_min.cdf(x, a, scale=b)
+        assert_allclose(c, -special.expm1(-0.25))
+
+        lc = stats.weibull_min.logcdf(x, a, scale=b)
+        assert_allclose(lc, np.log(-special.expm1(-0.25)))
+
+        s = stats.weibull_min.sf(x, a, scale=b)
+        assert_allclose(s, np.exp(-0.25))
+
+        ls = stats.weibull_min.logsf(x, a, scale=b)
+        assert_allclose(ls, -0.25)
+
+        # Also test using a large value x, for which computing the survival
+        # function using the CDF would result in 0.
+        s = stats.weibull_min.sf(30, 2, scale=3)
+        assert_allclose(s, np.exp(-100))
+
+        ls = stats.weibull_min.logsf(30, 2, scale=3)
+        assert_allclose(ls, -100)
+
+        # weibull_max
+        x = -1.5
+
+        p = stats.weibull_max.pdf(x, a, scale=b)
+        assert_allclose(p, np.exp(-0.25)/3)
+
+        lp = stats.weibull_max.logpdf(x, a, scale=b)
+        assert_allclose(lp, -0.25 - np.log(3))
+
+        c = stats.weibull_max.cdf(x, a, scale=b)
+        assert_allclose(c, np.exp(-0.25))
+
+        lc = stats.weibull_max.logcdf(x, a, scale=b)
+        assert_allclose(lc, -0.25)
+
+        s = stats.weibull_max.sf(x, a, scale=b)
+        assert_allclose(s, -special.expm1(-0.25))
+
+        ls = stats.weibull_max.logsf(x, a, scale=b)
+        assert_allclose(ls, np.log(-special.expm1(-0.25)))
+
+        # Also test using a value of x close to 0, for which computing the
+        # survival function using the CDF would result in 0.
+        s = stats.weibull_max.sf(-1e-9, 2, scale=3)
+        assert_allclose(s, -special.expm1(-1/9000000000000000000))
+
+        ls = stats.weibull_max.logsf(-1e-9, 2, scale=3)
+        assert_allclose(ls, np.log(-special.expm1(-1/9000000000000000000)))
+
 
 class TestRdist(TestCase):
     @dec.slow