sample_beta ill-conditioned for small values of alpha, beta about distributions HOT 4 CLOSED

This was discovered by @dpfau

from distributions.

A simple hack that's working for me is to add the line:

if (x==0 && y==0) return sample_bernoulli(rng, alpha / (alpha + beta)) ? 1.0 : 0.0;

right before the return statement. Johnk's algorithm would suffer from the same problem as the unmodified sampler here - if alpha = 0.001 and U = 0.5, U^(1/alpha) would underflow. You can see for yourself by trying out this gist: https://gist.github.com/dpfau/5ea001f82171ff8d9e99

from distributions.

Incidentally, numpy also produces NaNs quite often when alpha = beta = 0.001!

from distributions.

fixed by #90

from distributions.