Repository: pymc-devs/pymc
Pull request: Cap the values that Beta.random can generate.
Closes #3898.

To prevent beta random variates to be equal to 0 or 1, we cap the values to `np.nextafter(0,1, dtype)` and `np.nextafter(1,0, dtype)`, which give the next floating point nearby zero and one respectively, depending on the precision of the floating point representation.

Depending on what your PR does, here are a few things you might want to address in the description:
+ [X] what are the (breaking) changes that this PR makes?
Samples returned by `Beta.random` can never be equal to 0 or 1.

+ [x] important background, or details about the implementation

+ [X] are the changes—especially new features—covered by tests and docstrings?
Yes

+ [x] consider adding/updating relevant example notebooks
+ [X] right before it's ready to merge, mention the PR in the RELEASE-NOTES.md


Problem statement:
Difference in beta support between pymc3 and scipy results in error in sample_prior_predictive()
Description of your problem

**Please provide a minimal, self-contained, and reproducible example.**
```python
import numpy as np
import pymc3 as pm

idx = np.random.randint(0,1000,10000)

with pm.Model() as model:
    mu = pm.Beta('mu', alpha=.01, beta=.01, shape=1000)
    nu = pm.Gamma('nu', mu=30, sigma=5, shape=1000)
    theta = pm.Beta('theta', alpha=mu[idx] * nu[idx],  beta=(1-mu[idx])*nu[idx], shape=10000)
    
    
with model:
    prior_sample = pm.sample_prior_predictive()
```

**Please provide the full traceback.**
```python
Traceback (most recent call last):

  File "<ipython-input-39-03487933ffab>", line 14, in <module>
    prior_sample = pm.sample_prior_predictive()

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/sampling.py", line 1495, in sample_prior_predictive
    values = draw_values([model[name] for name in names], size=samples)

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/distributions/distribution.py", line 617, in draw_values
    value = _draw_value(next_,

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/distributions/distribution.py", line 787, in _draw_value
    return param.random(point=point, size=size)

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/model.py", line 56, in __call__
    return getattr(self.obj, self.method_name)(*args, **kwargs)

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/distributions/continuous.py", line 1264, in random
    return generate_samples(stats.beta.rvs, alpha, beta,

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/pymc3/distributions/distribution.py", line 960, in generate_samples
    samples = generator(size=dist_bcast_shape, *args, **kwargs)

  File "/home/zack/miniconda3/envs/_/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py", line 966, in rvs
    raise ValueError("Domain error in arguments.")

ValueError: Domain error in arguments.
```

**Please provide any additional information below.**

I am using beta as a hyperprior for other beta distributions and parameterizing based on mean and sample size. As a result, if the random variable from the hyper prior is 0 or 1, alpha or beta will end up 0, which is outside their domain. pymc3.Beta has support on (0,1), so I haven't seen an issue with sampling posteriors. When running sample_prior_predictive(), I get a domain error, which I believe is because scipy.stats.beta().rvs() will happily return 0 or 1.

Maybe this is just because of truncation in float64 (1 is far more common than 0), but I noticed the scipy.stats.beta documentation does state support on [0,1]. Not sure whether or not that makes it a scipy issue instead, but thought it was worth mentioning here.

## Versions and main components

* PyMC3 Version: 3.8
* Theano Version: 1.0.4
* Python Version: 3.8.1
* Operating system: Ubuntu 16.04
* How did you install PyMC3: (conda/pip) conda-forge

Hints from issue discussion:
Hi Zachary,
Thanks for reporting! This seems like a bug to me, as I'm under the impression that the Beta distribution's support should be [0, 1] and alpha and beta are both >= 0 (see [here](https://betanalpha.github.io/assets/case_studies/probability_densities.html#24_the_beta_family) for instance). But maybe someone knows why alpha and beta must be > 0 in `pm.Beta`?
Alpha and beta have to be > 0. The pdf of Beta contains &Gamma;(&alpha;) and  &Gamma;(&beta;). The gamma function is infinite at zero.
Some parameterizations/transformations have support on the closed interval, but for ours it should be open.
> Some parameterizations/transformations have support on the closed interval, but for ours it should be open.

@fonnesbeck Shouldn't the beta only be 0 when α = β = 1?  This is the definition in the doc page for `scipy.stats.beta`:

![image](https://user-images.githubusercontent.com/3274/81123528-e4ce9980-8ef8-11ea-84cd-1255fe8de9af.png)

And if x = 0, then x^{a - 1} can only be non-zero when a = 1.

@zack-kimble 's example shows (and I checked this) that `scipy.stats.beta(a, b).rvs()` can return 0 when `a` = `b` = 0.1, which seems wrong.

This is not sufficiently my thing that I don't feel confident reporting this to the scipy.stats project -- @zack-kimble would you do that?

I also don't know what we should do about this. Writing a high performance RNG for betas using numpy is really not my thing, but if someone has one that doesn't have this pathology, I'd be happy to try to plug it in.

Alternatively, should we just return `np.amax(samples, np.finfo(float).eps)` where we previously returned `samples`?  That's a quick and brutal fix.
When alpha and beta are 1, the beta distribution just becomes a uniform distribution.
@fonnesbeck I have replicated this issue, though, and we can get x = 0 when alpha and beta are 0.1, which is definitely wrong.

@zack-kimble We are using `scipy.stats.beta` so I think yes, you *should* report this as a bug to them.  I was easily able to replicate the issue invoking their beta random generator with a=b=0.01 as in your example, without any involvement of PyMC3.  Please keep us posted!
@zack-kimble, we've managed to identify that this is a problem caused by `numpy.random.beta`, and have opened an issue and a PR over at numpy. If the numpy devs agree with us that the beta's pdf domain should be the open (0, 1) interval, the fix will only become available once numpy makes a new release.
@lucianopaz, @rpgoldman  Thanks for tracking it back to numpy and opening a PR there. Sorry I didn't follow up myself there earlier!
After some discussion with the numpy devs, their random generator will remain as is, and we will clip the values of the samples it generates to lie in the open (0, 1) interval in pymc3. I hope to be able to open a PR tomorrow.

Commit message: Fix the PR link in the release notes

Patch to review:
diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index d593ad2141..6ebf36917c 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -25,6 +25,7 @@
 - End of sampling report now uses `arviz.InferenceData` internally and avoids storing
   pointwise log likelihood (see [#3883](https://github.com/pymc-devs/pymc3/pull/3883)).
 - The multiprocessing start method on MacOS is now set to "forkserver", to avoid crashes (see issue [#3849](https://github.com/pymc-devs/pymc3/issues/3849), solved by [#3919](https://github.com/pymc-devs/pymc3/pull/3919)).
+- Forced the `Beta` distribution's `random` method to generate samples that are in the open interval $(0, 1)$, i.e. no value can be equal to zero or equal to one (issue [#3898](https://github.com/pymc-devs/pymc3/issues/3898) fixed by [#3924](https://github.com/pymc-devs/pymc3/pull/3924)).
 
 ### Deprecations
 - Remove `sample_ppc` and `sample_ppc_w` that were deprecated in 3.6.
diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py
index c1c511e835..dfa82dc70b 100644
--- a/pymc3/distributions/continuous.py
+++ b/pymc3/distributions/continuous.py
@@ -32,6 +32,7 @@
 from .dist_math import (
     alltrue_elemwise, betaln, bound, gammaln, i0e, incomplete_beta, logpow,
     normal_lccdf, normal_lcdf, SplineWrapper, std_cdf, zvalue,
+    clipped_beta_rvs,
 )
 from .distribution import (Continuous, draw_values, generate_samples)
 
@@ -1290,7 +1291,7 @@ def random(self, point=None, size=None):
         """
         alpha, beta = draw_values([self.alpha, self.beta],
                                   point=point, size=size)
-        return generate_samples(stats.beta.rvs, alpha, beta,
+        return generate_samples(clipped_beta_rvs, alpha, beta,
                                 dist_shape=self.shape,
                                 size=size)
 
diff --git a/pymc3/distributions/dist_math.py b/pymc3/distributions/dist_math.py
index 81c595c921..6c2c772453 100644
--- a/pymc3/distributions/dist_math.py
+++ b/pymc3/distributions/dist_math.py
@@ -33,6 +33,11 @@
 
 f = floatX
 c = - .5 * np.log(2. * np.pi)
+_beta_clip_values = {
+    dtype: (np.nextafter(0, 1, dtype=dtype), np.nextafter(1, 0, dtype=dtype))
+    for dtype in ["float16", "float32", "float64", "float128"]
+}
+
 
 
 def bound(logp, *conditions, **kwargs):
@@ -548,3 +553,45 @@ def incomplete_beta(a, b, value):
             tt.and_(tt.le(b * value, one), tt.le(value, 0.95)),
             ps,
             t))
+
+
+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
+    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).
+    b : float or array_like of floats
+        Beta, 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),
+        a single value is returned if ``a`` and ``b`` are both scalars.
+        Otherwise, ``np.broadcast(a, b).size`` samples are drawn.
+    dtype : str or dtype instance
+        The floating point precision that the samples should have. This also
+        determines the value that will be used to shift any samples returned
+        by the numpy random number generator that are zero or one.
+    
+    Returns
+    -------
+    out : ndarray or scalar
+        Drawn samples from the parameterized beta distribution. The numpy
+        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)``.
+
+    """
+    out = np.random.beta(a, b, size=size).astype(dtype)
+    lower, upper = _beta_clip_values[dtype]
+    out[out == 0] = lower
+    out[out == 1] = upper
+    return out
\ No newline at end of file
diff --git a/pymc3/tests/test_dist_math.py b/pymc3/tests/test_dist_math.py
index b55e9792bc..f54f91bc2e 100644
--- a/pymc3/tests/test_dist_math.py
+++ b/pymc3/tests/test_dist_math.py
@@ -24,7 +24,9 @@
 from ..theanof import floatX
 from ..distributions import Discrete
 from ..distributions.dist_math import (
-    bound, factln, alltrue_scalar, MvNormalLogp, SplineWrapper, i0e)
+    bound, factln, alltrue_scalar, MvNormalLogp, SplineWrapper, i0e,
+    clipped_beta_rvs,
+)
 
 
 def test_bound():
@@ -216,3 +218,11 @@ def test_grad(self):
         utt.verify_grad(i0e, [-2.])
         utt.verify_grad(i0e, [[0.5, -2.]])
         utt.verify_grad(i0e, [[[0.5, -2.]]])
+
+
+@pytest.mark.parametrize("dtype", ["float16", "float32", "float64", "float128"])
+def test_clipped_beta_rvs(dtype):
+    # Verify that the samples drawn from the beta distribution are never
+    # equal to zero or one (issue #3898)
+    values = clipped_beta_rvs(0.01, 0.01, size=1000000, dtype=dtype)
+    assert not (np.any(values == 0) or np.any(values == 1))
\ No newline at end of file


Code review comments:
diff_hunk: @@ -548,3 +553,45 @@ def incomplete_beta(a, b, value):
             tt.and_(tt.le(b * value, one), tt.le(value, 0.95)),
             ps,
             t))
+
+
+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
+    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).
+    b : float or array_like of floats
+        Beta, 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),
+        a single value is returned if ``a`` and ``b`` are both scalars.
+        Otherwise, ``np.broadcast(a, b).size`` samples are drawn.
+    dtype : str or dtype instance
+        The floating point precision that the samples should have. This also
+        determines the value that will be used to shift any samples returned
+        by the numpy random number generator that are zero or one.
+    
+    Returns
+    -------
+    out : ndarray or scalar
+        Drawn samples from the parameterized beta distribution. The numpy
+        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)``.
path: pymc3/distributions/dist_math.py
text: @user1:
Typo: "... equal to 1 ~~if~~ **is** shifted to ..."

diff_hunk: @@ -548,3 +553,45 @@ def incomplete_beta(a, b, value):
             tt.and_(tt.le(b * value, one), tt.le(value, 0.95)),
             ps,
             t))
+
+
+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
+    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).
+    b : float or array_like of floats
+        Beta, positive (>0).
path: pymc3/distributions/dist_math.py
text: @user1:
Nitpicking: _strictly_ positive