### Java libraries for (non-uniform) random number simulation

Anyone writing serious Monte Carlo (and MCMC) codes relies on having a very good and fast (uniform) random number generator and associated functions for generation of non-uniform random quantities, such as Gaussian, Poisson, Gamma, etc. In a previous post I showed how to write a simple Gibbs sampler in four different languages. In C (and C++) random number generation is easy for most scientists, as the (excellent) GNU Scientific Library (GSL) provides exactly what most people need. But it wasn’t always that way… I remember the days before the GSL, when it was necessary to hunt around on the net for bits of C code to implement different algorithms. Worse, it was often necessary to hunt around for a bit of free FORTRAN code, and compile that with an F77 compiler and figure out how to call it from C. Even in the early Alpha days of the GSL, coverage was patchy, and the API changed often. Bad old days… But those days are long gone, and C programmers no longer have to worry about the problem of random variate generation – they can safely concentrate on developing their interesting new algorithm, and leave the rest to the GSL. Unfortunately for Java programmers, there isn’t yet anything quite comparable to the GSL in Java world.

I pretty much ignored Java until Java 5. Before then, the language was too limited, and the compilers and JVMs were too primitive to really take seriously for numerical work. But since the launch of Java 5 I’ve been starting to pay more interest. The language is now a perfectly reasonable O-O language, and the compilers and JVMs are pretty good. On a lot of benchmarks, Java is really quite comparable to C/C++, and Java is nicer to code, and has a lot of impressive associated technology. So if there was a math library comparable to the GSL, I’d be quite tempted to jump ship to the Java world and start writing all of my Monte Carlo codes in Java. But there isn’t. At least not yet.

When I first started to take Java seriously, the only good math library with good support for non-uniform random number generation was COLT. COLT was, and still is, pretty good. The code is generally well-written, and fast, and the documentation for it is reasonable. However, the structure of the library is very idiosyncratic, the coverage is a bit patchy, and there doesn’t ever seem to have been a proper development community behind it. It seems very much to have been a one-man project, which has long since stagnated. Unsurprisingly then, COLT has been forked. There is now a Parallel COLT project. This project is continuing the development of COLT, adding new features that were missing from COLT, and, as the name suggests, adding concurrency support. Parallel COLT is also good, and is the main library I currently use for random number generation in Java. However, it has obviously inherited all of the idiosyncrasies that COLT had, and still doesn’t seem to have a large and active development community associated with it. There is no doubt that it is an incredibly useful software library, but it still doesn’t really compare to the GSL.

I have watched the emergence of the Apache Commons Math project with great interest (not to be confused with Uncommons Math – another one-man project). I think this project probably has the greatest potential for providing the Java community with their own GSL equivalent. The Commons project has a lot of momentum, the Commons Math project seems to have an active development community, and the structure of the library is more intuitive than that of (Parallel) COLT. However, it is early days, and the library still has patchy coverage and is a bit rough around the edges. It reminds me a lot of the GSL back in its Alpha days. I’d not bothered to even download it until recently, as the random number generation component didn’t include the generation of gamma random quantities – an absolutely essential requirement for me. However, I noticed recently that the latest release (2.2) did include gamma generation, so I decided to download it and try it out. It works, but the generation of gamma random quantities is very slow (around 50 times slower than Parallel COLT). This isn’t a fundamental design flaw of the whole library – generating Gaussian random quantities is quite comparable with other libraries. It’s just that an inversion method has been used for gamma generation. All efficient gamma generators use a neat rejection scheme. In case anyone would like to investigate for themselves, here is a complete program for gamma generation designed to be linked against Parallel COLT:

import java.util.*; import cern.jet.random.tdouble.*; import cern.jet.random.tdouble.engine.*; class GammaPC { public static void main(String[] arg) { DoubleRandomEngine rngEngine=new DoubleMersenneTwister(); Gamma rngG=new Gamma(1.0,1.0,rngEngine); long N=10000; double x=0.0; for (int i=0;i<N;i++) { for (int j=0;j<1000;j++) { x=rngG.nextDouble(3.0,25.0); } System.out.println(x); } } }

and here is a complete program designed to be linked against Commons Math:

import java.util.*; import org.apache.commons.math.*; import org.apache.commons.math.random.*; class GammaACM { public static void main(String[] arg) throws MathException { RandomDataImpl rng=new RandomDataImpl(); long N=10000; double x=0.0; for (int i=0;i<N;i++) { for (int j=0;j<1000;j++) { x=rng.nextGamma(3.0,1.0/25.0); } System.out.println(x); } } }

The two codes do the same thing (note that they parameterise the gamma distribution differently). Both programs work (they generate variates from the same, correct, distribution), and the Commons Math interface is slightly nicer, but the code is *much* slower to execute. I’m still optimistic that Commons Math will one day be Java’s GSL, but I’m not giving up on Parallel COLT (or C, for that matter!) just yet…

Maybe you can consider the suite of random number generators from SuanShu:

http://www.numericalmethod.com/javadoc/suanshu/com/numericalmethod/suanshu/stats/random/RNG.html

It isn’t free software, and doesn’t include a gamma generator anyway.

You can create any generator you like using the Inverse Sampling method:

http://www.numericalmethod.com/javadoc/suanshu/com/numericalmethod/suanshu/stats/random/distribution/InverseTransformSampling.html

They already have the Gamma distribution provided:

http://www.numericalmethod.com/javadoc/suanshu/com/numericalmethod/suanshu/stats/distribution/univariate/GammaDistribution.html

It is probably free for your purpose: academic use.

As explained in the post, the reason the sampler (currently) in Commons Math is so slow is that it uses an inversion method. That is just not a good way to simulate gamma random quantities. Much more efficient algorithms exist, such as implemented in Parallel COLT.

Actually, SuanShu has a suite of gamma random number generators for you to choose from. Hope you like it!

http://www.numericalmethod.com/javadoc/suanshu/com/numericalmethod/suanshu/stats/random/distribution/gamma/package-summary.html

Hi, Darren.

I am developing my Randomness Framework.

http://code.google.com/p/randomness-framework/

This is collection of java RNG with propeties make it suitable for server-side applications. When RF will be completed, i am plan to extend it into simulation server and include varyous distribution’s.

SSJ is another option:

http://www.iro.umontreal.ca/~simardr/ssj/indexe.html

In my tests it was slower than Colt, but it is an interesting library with many useful features.

Hi Darren

Following our discussion last week (very nice, thanks again!), I’ve fiddled with Java, and it indeed has come a looong way from my last try on Java 1.2!

Regarding the above, I’ve tried your code above with Apache Commons Math 3 (same API, just change the import), which now uses a rejection algorithm for their Gamma., and while the source shows it’s as elaborate as PC’s algorithm its runtime is nothing like ACM2 . On my laptop:

ACM2: 107 seconds (!)

ACM3: Between 2.4 seconds and 6 seconds — varies quite largely.

PC: 1.8 seconds.

Yes, indeed. When I wrote this post, I also filed this as an issue with the ACM team, and they fixed it in ACM3 within a month – quite impressive, I think!

I recently found jMEF at:

http://www.lix.polytechnique.fr/~nielsen/MEF/

which implements mixtures of exponential families and has an MIT license.

The information on jMEF lists these distributions: “Gaussian (generic, isotropic Gaussian, diagonal Gaussian, rectified Gaussian or Wald distributions, lognormal), Poisson, Bernoulli, binomial, multinomial, Laplacian, Gamma (incl. chi-squared), Beta, exponential, Wishart, Dirichlet, Rayleigh, probability simplex, negative binomial distribution, Weibull, von Mises, Pareto distributions, skew logistic, etc.”

I do not know how it compares to ACM3 or PCOLT, but the src code looks clean/elegant.

I’m the author of jMEF, thank you for your compliment 🙂 I’ve spent a lot of time trying to make it clean. Mixture models are not that easy to code, especially when you try to provide a generic approach that allow you to estimate any kind of mixture of exponential family.

The link to jMEF changed. It’s now on GitHub: http://vincentfpgarcia.github.com/jMEF/

I believe a year after this post, Apache Commons Math no longer uses inversion sampling:

https://issues.apache.org/jira/browse/MATH-774

http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/random/RandomDataGenerator.html#nextGamma%28double,%20double%29