Generating Stable Distributions

The following are some C-code fragments you may find useful in generating values from stable distributions. These all make use of a routine called "randomfl()", which should be a small-space random number generator which returns a double distributed uniformly in the range 0.0 to 1.0. Use your favourite random number generator for this function. The following routines are useful depending on whether you want to compute the L1 distance, L2, Lp for arbitrary p between 0 and 2, or "L0" (Lp with p close to zero).

Code written by Graham Cormode except where otherwise acknowledged. This is made freely available for all to use, but please notify me if you find it of use, I would like to know where it is being used. Email graham a dimacs.rutgers.edu.
Processing Massive Data Sets home page

Cauchy distribution, p-stable with p=1

long double cauchy() {

  // return a value from the cauchy distribution
  // this is distributed Cauchy(1,0)

  return(tan(PI*(randomfl() - 0.5)));
}

Gaussian Distribution, p-stable with p=2

long double gasdev() {

  // Pick random values distributed N(0,1) using the Box-Muller transform
  // Taken from Numerical Recipes in C 
  // picks two at a time, returns one per call (buffers the other)

  static int iset=0;
  static long double gset;
  long double fac,rsq,v1,v2;
  if (iset == 0) {
    do {
      v1 = 2.0*randomfl()-1.0;
      v2 = 2.0*randomfl()-1.0;
      rsq=v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    
    fac = sqrt((long double) -2.0*log((long double)rsq)/rsq);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  }
  else {
    iset = 0;
    return gset;
  }
}

Stable distributions for arbitrary p

Call this with parameter p to generate values from a p-stable distribution.

double stabledevd(float alpha) {

  // we set beta = 0 by analogy with the normal and cauchy case
  // in order to get symmetric stable distributions. 

  double theta, W, holder, left, right;

  theta=PI*(randomfl() - 0.5);
  W = -log(randomfl()); // takes natural log

  // if beta == 0 then c(alpha,beta)=1; theta_0 = 0
  // expression reduces to sin alpha.theta / (cos theta) ^1/alpha
  //  * (cos (theta - alpha theta)/W) ^(1-alpha)/alpha

  left = (sin(alpha*theta)/pow(cos(theta), 1.0/alpha));
  right= pow(cos(theta*(1.0 - alpha))/W, ((1.0-alpha)/alpha));
  holder=left*right;
  return(holder);
}

Faster code for p close to zero

The following code replaces the stable distribution with a distribution that is in the domain of attraction of a stable distribution with parameter alpha. Thanks to John Nolan for suggesting this approach and some helpful comments. See "Comparing Data Streams Using Hamming Norms" in IEEE Transactions on Knowledge and Data Engineering, 2003.

double altstab(double alpha) 
{
  double u,v,result;
  
  u=randomfl();
  v=randomfl();
  result=pow(u,-1.0/alpha);
  if (v<0.5) result=-result;
  
  return(result);
}

(Commercial use prohibited without author's permission)