MAINT: enable rdist._cdf, making it a lot faster. Closes scipygh-1285.

Issue for hyp2f1 that is worked around is discussed in scipygh-1561.

scipy/stats/distributions.py

@@ -1,5 +1,3 @@
-# Functions to implement several important functions for
-#   various Continous and Discrete Probability Distributions
 #
 # Author:  Travis Oliphant  2002-2011 with contributions from
 #          SciPy Developers 2004-2011
@@ -5114,7 +5112,6 @@ def _ppf(self, q, c):
 powernorm = powernorm_gen(name='powernorm', shapes="c")
 
 
-# FIXME: PPF does not work.
 class rdist_gen(rv_continuous):
     """An R-distributed continuous random variable.
 
@@ -5132,15 +5129,21 @@ class rdist_gen(rv_continuous):
 
     """
     def _pdf(self, x, c):
-        return np.power((1.0-x*x),c/2.0-1) / special.beta(0.5,c/2.0)
+        return np.power((1.0 - x**2), c / 2.0 - 1) / special.beta(0.5, c / 2.0)
+
+    def _cdf(self, x, c):
+        term1 = x / special.beta(0.5, c / 2.0)
+        res = 0.5 + term1 * special.hyp2f1(0.5, 1 - c / 2.0, 1.5, x**2)
+        # There's an issue with hyp2f1, it returns nans near x = +-1, c > 100.
+        # Use the generic implementation in that case.  See gh-1285 for
+        # background.
+        if any(np.isnan(res)):
+            return rv_continuous._cdf(self, x, c)
 
-    def _cdf_skip(self, x, c):
-        # error in special.hyp2f1 for some values see tickets 758, 759
-        return 0.5 + x/special.beta(0.5,c/2.0) * \
-               special.hyp2f1(0.5,1.0-c/2.0,1.5,x*x)
+        return res
 
     def _munp(self, n, c):
-        return (1-(n % 2))*special.beta((n+1.0)/2,c/2.0)
+        return (1 - (n % 2)) * special.beta((n + 1.0) / 2, c / 2.0)
 rdist = rdist_gen(a=-1.0, b=1.0, name="rdist", shapes="c")
 
 

scipy/stats/tests/test_continuous_basic.py

@@ -102,8 +102,6 @@
     ['powernorm', (4.4453652254590779,)],
     ['rayleigh', ()],
     ['rdist', (0.9,)],   # feels also slow
-#    ['rdist', (3.8266985793976525,)],  #veryslow, especially rvs
-    #['rdist', (541.0,)],   # from ticket #758    #veryslow
     ['recipinvgauss', (0.63004267809369119,)],
     ['reciprocal', (0.0062309367010521255, 1.0062309367010522)],
     ['rice', (0.7749725210111873,)],

scipy/stats/tests/test_continuous_extra.py

@@ -123,6 +123,15 @@ def test_erlang_runtimewarning():
         npt.assert_allclose(result_erlang, result_gamma, rtol=1e-3)
 
 
+@npt.dec.slow
+def test_rdist_cdf_gh1285():
+    # check workaround in rdist._cdf for issue gh-1285.
+    distfn = stats.rdist
+    values = [0.001, 0.5, 0.999]
+    npt.assert_almost_equal(distfn.cdf(distfn.ppf(values, 541.0), 541.0),
+                            values, decimal=5)
+
+
 if __name__ == "__main__":
     import nose
     #nose.run(argv=['', __file__])