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 1


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).


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;

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.