# Random number generation for gaussian, cauchy and levy

I am conducting my thesis on Evolutionary Algorithm. To perform various type of mutation I have to generate random numbers from Gaussian, Levy, Cauchy distributions. How can I implement this? I need pseudocode or algorithm. Thanks in advance.

## 1 Answer

First of all, do you really need to implement these RNGs at a low level of code? If you are going to be, e.g., writing C code, you can link to the GNU Scientific Library, which has all three of those distributions and a lot more besides. Available on both Linux and Windows, and Mac too I'm sure! Higher level languages like R also have RNGs built in: rnorm for the Gaussian, rlevy (in the rmutil package) for the Levy, and rcauchy for the Cauchy in R.

To actually answer your question, assuming you have a good uniform RNG, and know how to generate exponential variates, I've attached the GSL C source for the Levy, which includes the Cauchy and Gaussian as special cases, and some extra comments in the code (although if speed is an issue, you can do better for the Cauchy and Gaussian.) It's not 100% transparent, as it calls other GSL routines, but it should be obvious what's going on with those routines.

/* The stable Levy probability distributions have the form

p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)

with 0 < alpha <= 2.

For alpha = 1, we get the Cauchy distribution
For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.

Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
Simulation". The original reference given there is,

J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
simulating stable random variates". Journal of the American
Statistical Association, JASA 71 340-344 (1976).

*/

double
gsl_ran_levy (const gsl_rng * r, const double c, const double alpha)
{
double u, v, t, s;

/* gsl_rng_uniform_pos(r) returns a U(0,1) variate */

u = M_PI * (gsl_rng_uniform_pos (r) - 0.5);

if (alpha == 1)               /* cauchy case */
{
t = tan (u);
return c * t;
}

do
{
v = gsl_ran_exponential (r, 1.0); /* Just use -log(u) where u~U(0,1) */
}
while (v == 0);

if (alpha == 2)               /* gaussian case */
{
t = 2 * sin (u) * sqrt(v);
return c * t;
}

/* general case */

t = sin (alpha * u) / pow (cos (u), 1 / alpha);
s = pow (cos ((1 - alpha) * u) / v, (1 - alpha) / alpha);

return c * t * s;
}

• Could you please post the pseudocode! GSL documentation is very difficult to understand
– user
Jan 1, 2012 at 7:06
• Here's the java version of my code pastebin.com/RcMhFsuZ. Please have a look. I would really appreciate it
– user
Jan 1, 2012 at 8:08
• @crucifiedsoul - your code looks fine to me. Next time I'll give you Java code! Jan 1, 2012 at 16:11
• According to your code what will be the mean, standard deviation of Gaussian Distribution?
– user
Jan 8, 2012 at 8:35
• The mean will be zero, the standard deviation will be c*sqrt(2). Jan 8, 2012 at 15:54