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

Hi,

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?

http://forums.nvidia.com/index.php?showtopic=150988&pid=956395&mode=threaded&start=0#entry956395
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.

                            Bill

    Dr. W. B. Langdon,
    Department of Computer Science,
    University College London
    Gower Street, London WC1E 6BT, UK
    http://www.cs.ucl.ac.uk/staff/W.Langdon/

CIGPU 2011 CIGPU-2012 WCCI-2012 IJCNN-2012, CEC2012
A Field Guide to Genetic Programming
http://www.gp-field-guide.org.uk/
GECCO 2011 GECCO 2011
GP EM Genetic Programming and Evolvable Machines | Home
GP Bibliography The Genetic Programming 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:

http://developer.download.nvidia.com/compute/cuda/3_2/toolkit/docs/CURAND_Library.pdf

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){
    do{
        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) 
      break;

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

}
}