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