Currently the `BetaDistribution#sample` uses inverse CDF method, which is quite slow for sampling-intensive computations. I've implemented a method from the R. C. H. Cheng paper and it seems to work much better. Here's a simple microbenchmark:
Should I submit a patch?
R. C. H. Cheng (1978). Generating beta variates with nonintegral shape parameters. Communications of the ACM, 21, 317–322.