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?
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/
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?
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;