Beta distribution pdf returns inf using output of Beta rvs

[gobbedy]

Hello,

The beta sampler (beta.rvs) sometimes returns a 1.0 value.

In turn when passed into the beta pdf function (beta.pdf) this 1.0 return an inf values.

Does the 1.0 show up because of some float64 rounding error?

And if so, can I work around this somehow? I need to pass the output of rvs to the pdf as part of some math operations in my code.

The slight over-representation of 1.0 caused by float64 precision rounding (assuming that is the case) is not a problem, but the inf causes many simulations to fail.

Does the beta distribution support long double?

I would be grateful for any help.

>>> rvs=beta.rvs(1000,0.1,size=20)
>>> rvs
array([1.                , 0.9999999999485164, 0.9999997435379178,
       0.9999999811610324, 0.9999404607880107, 0.9999999895713109,
       0.9999980001104486, 0.9999999999999962, 0.9999999994453201,
       0.9999999999886771, 0.9999851060013317, 0.9997969567848644,
       0.9992273211069278, 0.9999999950528663, 0.9999999998670334,
       0.9996860592299935, 0.9998663752734042, 0.9987631795220868,
       0.9999999993056435, 0.9998308307442249])
>>> beta.pdf(rvs,1000,0.1)
array([                   inf, 3.8118673391078585e+08,
       1.7922874503637722e+05, 1.8796639719161065e+06,
       1.2545400875895550e+03, 3.2005786647235840e+06,
       2.8175179827495387e+04, 2.0065017738705874e+12,
       4.4874681965535089e+07, 1.4896231077731190e+09,
       4.5652823162722752e+03, 3.6033698030260530e+02,
       6.1245818869958853e+01, 6.2620906577958884e+06,
       1.6228235603163967e+08, 2.1789589412901623e+02,
       5.6281195096702629e+02, 2.5212998507201423e+01,
       3.6661923420780651e+07, 4.3928712163659350e+02])

[pvanmulbregt]

The PDF for the beta distribution is given by

PDF(x; alpha, beta) = x**(alpha-1) * (1-x)**(beta-1) / B(alpha, beta)

In the example above, alpha=1000, beta=0.1, so that as x approaches 1, the PDF is approximately (1-x)**(-0.9), which does indeed approach infinity.
The implementation of beta._rvs(self, a, b) calls self._random_state.beta(a, b, self._size) which appears to be a numpy function to generate the random sample.

[IandJMSmith]

If you use
rvs=beta.rvs(0.1,1000,size=20)
then
beta.pdf(rvs,0.1,1000)
will give more accurate versions of the values you would expect from your your original example.
In particular you would be very unlucky to get an inf value.

Basically, you are now generating numbers very close to zero instead of numbers very close to 1. This avoids the (1-x) calculation which due to rounding errors currently delivers 0 rather than someting close to 0.

[mdhaber]

It appears from the definition that the PDF at x = 1 is indeed inf. (R agrees, if it's any help. dbeta(0, 0.1, 1000)). If I understand correctly, the post was not about the rvs, but if it were, it should probably be reported to NumPy; np.random.beta(1000, 0.1, size=20) gives the same samples.
I don't see a bug, so I'm closing. If I've misunderstood, please reopen @gobbedy. Thanks!