Pseudorandom Number Generator - Exponential Distri

Question

I would like to generate some pseudorandom numbers and up until now I've been very content with the .Net library's Random.Next(int min, int max) function. PRNGs of this variety are supposed to be using a Uniform distribution, but I would very much like to generate some numbers using an Exponential Distribution.

I'm programming in C#, although I'll accept pseudocode or C++, Java or the like.

Any suggestions / code snippets / algorithms / thoughts?

Answer 1

One interesting property of the exponential distribution: Consider an arrival process with exponential interarrival times. Take any period of time (t1, t2) and the arrivals in that period. Those arrivals are UNIFORMLY distributed between t1 and t2. (Sheldon Ross, Stochastic Processes).

If I have a pseudo-random number generator and, for some reason (e.g. my software can't calculate logs) you don't want to do the above transformation, but want an exponential r.v. with mean of 1.0.

You can :

1) Create 1001 U(0,1) random variables.

2) Sort in order

3) Subtract the second from the first, third from the second,... to get 1000 differences.

4) Those differences are exponential RVs with from a distribution with mean = 1.0.

Less efficient, I think, but a means to the same end.

Answer 2

This was what I used when faced with similar requirements:

// sorry.. pseudocode, mine was in Tcl:

int weighted_random (int max) {
    float random_number = rand();
    return floor(max - ceil( max * random_number * random_number))
}

Of course this is the formula of squaring the random number so you're generating a random number along a quadratic curve.

Answer 3

Since you have access to a uniform random number generator, generating a random number distributed with other distribution whose CDF you know is easy using the inversion method.

So, generate a uniform random number, u, in [0,1), then calculate x by:

x = log(1-u)/(−λ),

where λ is the rate parameter of the exponential distribution. Now, x is a random number with an exponential distribution. Note that log above is ln, the natural logarithm.

Answer 4

The Fundamental Theorem of Sampling holds that if you can normalize, integrate and invert the desired distribution you are home free.

If you have a desired distribution F(x) normalized on [a,b]. You compute

C(y) = \int_a^y F(x) dx

invert that to get C^{-1}, throw z uniformly on [0,1) and find

x_i = C^{-1}(z_i)

which will have the desired distribution.

---

In your case: F(x) = ke^{-kx} and I will assume that you want [0,infinity]. We get :

C(y) = 1 - e^{-ky}

which is invertable to give

x = -1/k  ln(1 - z)

for z thrown uniformly on [0,1).

---

But, frankly, using a well debugged library is smarter unless you're doing this for your own edification.

Answer 5

If you want good random numbers, consider linking to the gsl routines: http://www.gnu.org/software/gsl/. They have the routine gsl_ran_exponential. If you want to generate random numbers using a built-in generator with a uniform distribution on [0, 1) (e.g. u=Random.Next(0, N-1)/N, for some large N), then just use:

-mu * log (1-u)

See randist/exponential.c in the gsl source.

EDIT: just for comparison with some later answers - this is equivalent with mu = 1/lambda. mu here is the mean of the distribution, also called the scale parameter on the wikipedia page the OP linked to, and lambda is the rate parameter.

Answer 6

The open-source Uncommons Maths library by Dan Dyer provides random number generators, probability distributions, combinatorics and statistics for Java.

Among other valuable classes, ExponentialGenerator has essentially implemented the idea explained by @Alok Singhal. In its tutorial blog, a code snippet is given to simulate some random event that happened on average 10 times a minute:

final long oneMinute = 60000;
Random rng = new MersenneTwisterRNG();

// Generate events at an average rate of 10 per minute.
ExponentialGenerator gen = new ExponentialGenerator(10, rng);
boolean running = true;
while (true)
{
    long interval = Math.round(gen.nextValue() * oneMinute);
    Thread.sleep(interval);

    // Fire event here.
}

Of course, if you prefer to the time unit per second (instead of a minute here), you just need to set final long oneMinute = 1000.

Going deeper into the source code of the method nextValue() of ExponentialGenerator, you will find the so-called inverse transform sampling described in Generating_exponential_variates [wiki]:

public Double nextValue()
{
    double u;
    do
    {
        // Get a uniformly-distributed random double between
        // zero (inclusive) and 1 (exclusive)
        u = rng.nextDouble();
    } while (u == 0d); // Reject zero, u must be positive for this to work.
    return (-Math.log(u)) / rate.nextValue();
}  

P.S.: Recently I am using the Uncommons Maths library. Thanks Dan Dyer.