FIX: use scipy.stats and change ref_rand to point to clipped_beta_rvs

pymc3/distributions/dist_math.py

@@ -19,6 +19,7 @@
 '''
 import numpy as np
 import scipy.linalg
+import scipy.stats
 import theano.tensor as tt
 import theano
 from theano.scalar import UnaryScalarOp, upgrade_to_float_no_complex
@@ -558,17 +559,17 @@ def incomplete_beta(a, b, value):
 def clipped_beta_rvs(a, b, size=None, dtype="float64"):
     """Draw beta distributed random samples in the open :math:`(0, 1)` interval.
 
-    The samples are generated with ``numpy.random.beta``, but any value that
+    The samples are generated with ``scipy.stats.beta.rvs``, but any value that
     is equal to 0 or 1 will be shifted towards the next floating point in the
     interval :math:`[0, 1]`, depending on the floating point precision that is
     given by ``dtype``.
 
     Parameters
     ----------
     a : float or array_like of floats
-        Alpha, positive (>0).
+        Alpha, strictly positive (>0).
     b : float or array_like of floats
-        Beta, positive (>0).
+        Beta, strictly positive (>0).
     size : int or tuple of ints, optional
         Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
         ``m * n * k`` samples are drawn.  If size is ``None`` (default),
@@ -582,15 +583,15 @@ def clipped_beta_rvs(a, b, size=None, dtype="float64"):
     Returns
     -------
     out : ndarray or scalar
-        Drawn samples from the parameterized beta distribution. The numpy
+        Drawn samples from the parameterized beta distribution. The scipy
         implementation can yield values that are equal to zero or one. We
         assume the support of the Beta distribution to be in the open interval
         :math:`(0, 1)`, so we shift any sample that is equal to 0 to
         ``np.nextafter(0, 1, dtype=dtype)`` and any sample that is equal to 1
-        if shifted to ``np.nextafter(1, 0, dtype=dtype)``.
+        is shifted to ``np.nextafter(1, 0, dtype=dtype)``.
 
     """
-    out = np.random.beta(a, b, size=size).astype(dtype)
+    out = scipy.stats.beta.rvs(a, b, size=size).astype(dtype)
     lower, upper = _beta_clip_values[dtype]
     out[out == 0] = lower
     out[out == 1] = upper

pymc3/tests/test_distributions_random.py

@@ -22,6 +22,7 @@
 import theano
 
 import pymc3 as pm
+from pymc3.distributions.dist_math import clipped_beta_rvs
 from pymc3.distributions.distribution import (draw_values,
                                               _DrawValuesContext,
                                               _DrawValuesContextBlocker)
@@ -548,7 +549,7 @@ def ref_rand(size, mu, lam, alpha):
 
     def test_beta(self):
         def ref_rand(size, alpha, beta):
-            return st.beta.rvs(a=alpha, b=beta, size=size)
+            return clipped_beta_rvs(a=alpha, b=beta, size=size)
         pymc3_random(pm.Beta, {'alpha': Rplus, 'beta': Rplus}, ref_rand=ref_rand)
 
     def test_exponential(self):