Generating random samples from beta, gamma and multinomial distributions in the kernel


I’m working on a code that requires running the equivalent of the gsl functions gsl_ran_multinomial, gsl_ran_beta and gsl_ran_gamma on the kernel. These functions generate random samples from the multinomial, beta and gamma distributions correspondingly.
Is there a simple way to do that?
says how to generate a uniformly distributed psuedo random number quickly in CUDA
Numerical recipies in C has a section on converting uniform PRNGs to a number of other distributions.
The C implemention of their algorithms is available online.


    Dr. W. B. Langdon,
    Department of Computer Science,
    University College London
    Gower Street, London WC1E 6BT, UK

CIGPU 2011
A Field Guide to Genetic Programming
GECCO 2011
GP Bibliography

I would also suggest you take a look at the CURAND library that comes with the CUDA toolkit. It is very easy to use, and the XORWOW algorithm is very fast. (The period is shorter than Mersenne Twister, but 2^192 is good enough for most applications.) The documentation is here:

Regardless, basically all generators of non-uniform random number distributions start with a uniform RNG. You can take a look at the GSL source code directly (or Numerical Recipes) and in many cases, it should be fairly straightforward to port their function over to CUDA, replacing calls to the uniform RNG with curand_uniform().

I have ported over the implementation of exponential and Poisson distributions from CERN’s ROOT library, but unfortunately none of the distributions you mention.

Hello, it is easy to generate uniform random number using CURAND library. But when I try to apply it to generator of Gamma or Poisson distribution, I don’t know how to take advantage of GPU. How can I call a curand_uniform kernel inside my kernel to generate random number from Gamma distribution?

Hi, have you solved this question?

This is what I use at the moment.
It is the GSL function modified to work with CuRAND.

device double ran_gamma (curandState localState, const double a, const double b){
/* assume a > 0 */

if (a < 1){
double u = curand_uniform_double(&localState);
return ran_gamma (localState, 1.0 + a, b) * pow (u, 1.0 / a);

double x, v, u;
double d = a - 1.0 / 3.0;
double c = (1.0 / 3.0) / sqrt (d);

while (1){
        x = curand_normal_double(&localState);
        v = 1.0 + c * x;
    } while (v <= 0);

    v = v * v * v;
    u = curand_uniform_double(&localState);

    if (u < 1 - 0.0331 * x * x * x * x) 

    if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
return b * d * v;