CODEWITHSUNDEEP
c#randomdistribution
Wayne
Wayne

Reputation: 2989

Stable distribution random numbers?

How to generate random numbers with a stable distribution in C#? The Random class has uniform distribution. Many other code on the internet show normal distribution. But we need stable distribution meaning infinite variance, a.k.a fat-tailed distribution.

The reason is for generating realistic stock prices. In the real world, huge variations in prices are far more likely than in normal distributions.

Does someone know the C# code to convert Random class output into stable distribution?

Edit: Hmmm. Exact distribution is less critical than ensuring it will randomly generate huge sigma like at least 20 sigma. We want to test a trading strategy for resilience in a true fat tailed distribution which is exactly how stock market prices behave.

I just read about ZipFian and Cauchy due to comments. Since I must pick, let's go with Cauchy distribution but I will also try ZipFian to compare.

Upvotes: 4

Views: 1899

Answers (4)

Peter O.
Peter O.

Reputation: 32878

The following C# code generates a random number following a stable distribution given the shape parameters alpha and beta. I release it to the public domain under Creative Commons Zero.

public static double StableDist(Random rand, double alpha, double beta){
    if(alpha<=0 || alpha>2 || beta<-1 || beta>1)
        throw new ArgumentException();
    var halfpi=Math.PI*0.5;
    var unif=NextDouble(rand);
    while(unif == 0.0)unif=NextDouble(rand);
    unif=(unif - 0.5) * Math.PI;
    // Cauchy special case
    if(alpha==1 && beta==0)
        return Math.Tan(unif);
    var expo=-Math.Log(1.0 - NextDouble(rand));
    var c=Math.Cos(unif);
    if(alpha == 1){
        var s=Math.Sin(unif);
        return 2.0*((unif*beta+halfpi)*s/c -
            beta * Math.Log(halfpi*expo*c/(
                unif*beta+halfpi)))/Math.PI;
    }
    var z=-Math.Tan(halfpi*alpha)*beta;
    var ug=unif+Math.Atan(-z)/alpha;
    var cpow=Math.Pow(c, -1.0 / alpha);
    return Math.Pow(1.0+z*z, 1.0 / (2*alpha))*
        (Math.Sin(alpha*ug)*cpow)*
        Math.Pow(Math.Cos(unif-alpha*ug)/expo, (1.0-alpha) / alpha);
}

private static double NextDouble(Random rand){
    // The default NextDouble implementation in .NET (see
    // https://github.com/dotnet/corert/blob/master/src/System.Private.CoreLib/shared/System/Random.cs)
    // is very problematic:
    // - It generates a random number 0 or greater and less than 2^31-1 in a 
    // way that very slightly biases 2^31-2.
    // - Then it divides that number by 2^31-1.
    // - The result is a number that uses roughly only 32 bits of pseudorandomness, 
    // even though `double` has 53 bits in its significand.
    // To alleviate some of these problems, this method generates a random 53-bit
    // random number and divides that by 2^53.  Although this doesn't fix the bias
    // mentioned above (for the default System.Random), this bias may be of
    // negligible importance for most purposes not involving security.
    long x=rand.Next(0,1<<30);
    x<<=23;
    x+=rand.Next(0,1<<23);
    return (double)x / (double)(1L<<53);
}

In addition, I set forth pseudocode for the stable distribution in a separate article.

Upvotes: 0

Eric Lippert
Eric Lippert

Reputation: 660169

In general, the method is:

  • Choose a stable, fat-tailed distribution. Say, the Cauchy distribution.

  • Look up the quantile function of the chosen distribution.

For the Cauchy distribution, that would be p --> peak + scale * tan( pi * (p - 0.5) ).

  • And now you have a method of transforming uniformly-distributed random numbers into Cauchy-distributed random numbers.

Make sense? See

http://en.wikipedia.org/wiki/Inverse_transform_sampling

for details.

Caveat: It has been a long, long time since I took statistics.

UPDATE:

I liked this question so much I just blogged it: see

http://ericlippert.com/2012/02/21/generating-random-non-uniform-data/

My article exploring a few interesting examples of Zipfian distributions is here:

http://blogs.msdn.com/b/ericlippert/archive/2010/12/07/10100227.aspx

Upvotes: 7

David Heffernan
David Heffernan

Reputation: 613063

I have before me James Gentle's Springer volume on this topic, Random Number Generation and Monte Carlo Methods, courtesy of my statistician wife. It discusses the stable family on page 105:

The stable family of distributions is a flexible family of generally heavy-tailed distributions. This family includes the normal distribution at one extreme value of one of the parameters and the Cauchy at the other extreme value. Chambers, Mallows, and Stuck (1976) give a method for generating deviates from stable distributions. (Watch for some errors in the constants in the auxiliary function D2, for evaluating (ex-1)/x.) Their method is used in the IMSL libraries. For a symmetric stable distribution, Devroye (1986) points out that a faster method can be developed by exploiting the relationship of the symmetric stable to the Fejer-de la Vallee Poissin distribution. Buckle (1995) shows how to simulate the parameters of a stable distribution, conditional on the data.

Generating deviates from the generic stable distribution is hard. If you need to do this then I would recommend a library such as IMSL. I do not advise you attempt this yourself.

However, if you are looking for a specific distribution in the stable family, e.g. Cauchy, then you can use the method described by Eric, known as the probability integral transform. So long as you can write down the inverse of the distribution function in closed form then you can use this approach.

Upvotes: 0

LBushkin
LBushkin

Reputation: 131726

If you're interested in using the Zipfian distribution (which is often used when modeling processes from the sciences or social domains), you would do something along the lines of:

  1. Select your k (skew) for the distribution
  2. Precompute the domain of the cumulative distribution (this is just an optimization)
  3. Generate random values for the distribution by finding the nearest value from the domain

Sample Code:

List<int> domain = Enumerable.Range(0,1000);  // generate your domain
double skew  = 0.37; // select a skew appropriate to your domain
double sigma = domain.Aggregate(0.0d, (z,x) => x + 1.0 / Math.Pow(z+1, skew));
List<double> cummDist = domain.Select( 
      x => domain.Aggregate(0.0d, (z,y) => z + 1.0/Math.Pow(y, skew) * sigma));

Now you can generate random values by selecting the closest value from within the domain:

Random rand = new Random();
double seek = rand.NextDouble();
int searchIndex = cummDist.BinarySearch(seek);
// return the index of the closest value from the distribution domain
return searchIndex < 0 ? (~searchIndex)-1 : searchIndex-1;

You can, of course, generalize this entire process by factoring out the logic that materializes the domain of the distribution from the process that maps and returns a value from that domain.

Upvotes: 1

Related Questions