BUG: The beta distribution generates random values equal to 0

[tomicapretto]

Describe the issue:

Someone reported a test failing in Bambi (see bambinos/bambi#888).

The test was checking whether random values generated from a beta likelihood were always in the (0, 1) interval.

I managed to reproduce the example in PyMC.

Reproduceable code example:

import numpy as np
import pymc as pm

mu = np.array(
    [
        7.07551617e-12,
        3.13190732e-11,
        1.37305053e-10,
        5.71933452e-08,
        3.47710424e-07,
        7.36626506e-06,
        1.69159968e-04,
        4.39823566e-03,
        3.43223695e-01,
        3.41184856e-01,
    ]
)

kappa = np.array(
    [
        2.42025905,
        2.40677827,
        2.40394316,
        2.36733373,
        2.36652374,
        2.3536546,
        2.35645956,
        2.34899749,
        2.37637218,
        2.35059481,
    ]
)

alpha = mu * kappa
beta = (1 - mu) * kappa

draws = pm.draw(pm.Beta.dist(alpha, beta))
print(draws)
print((draws == 0).sum())

Error message:

PyMC version information:

import pymc as pm
import pytensor as pt

print(pm.__version__)
print(pt.__version__)

5.20.1
2.28.1

Context for the issue:

No response

[ricardoV94]

We actually have a clipped_beta here, but we are not using it as the rv_op. But are we sure we want to use it?

PyTensor is just defaulting to what numpy does. I guess numba does the same, and not sure about JAX.

pymc/pymc/distributions/continuous.py

Lines 1046 to 1052 in a9c5a85

	class BetaClippedRV(BetaRV):
	@classmethod
	def rng_fn(cls, rng, alpha, beta, size) -> np.ndarray:
	return np.asarray(clipped_beta_rvs(alpha, beta, size=size, random_state=rng))
	beta = BetaClippedRV()

pymc/pymc/distributions/continuous.py

Line 1129 in a9c5a85

rv_op = pytensor.tensor.random.beta

Boundary situations are always tricky with float precision. And the support is ambiguous? Wikipedia includes both: https://en.wikipedia.org/wiki/Beta_distribution

Our logp also accepts 0 or 1

[tomicapretto]

But not only the probability at x = {0, 1} is 0, but also the density (edit: it can also be undefined when alpha or beta are < 1). That's why I think they should be excluded.

(btw I agree the support is ambiguous)

[ricardoV94]

But not only the probability at x = {0, 1} is 0

Not if alpha=beta=1 ;)

[ricardoV94]

Anyway I would say if numpy is happy with it, we don't need to deviate here. Strong disagreement?

[tomicapretto]

Nope. I think it's totally fine. I can adjust the test on the Bambi side. Thanks!