Pseudorandom Number Generator - Exponential Distribution

I would like to generate some pseudorandom numbers and up until now开发者_StackOverflow社区 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?

---

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.

---

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.

---

This is the formula I found on Wikipedia :

T = -Ln(u) / λ

We create a random number with a uniform distribution (u) in [0,1]and we get x :

Random R = new Random();

double u = R. NextDouble();

double x = -Math.Log(u)/(λ);

---

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.

---

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.

---

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.

---

There is another way to generate an exponential(rate) random variate, although it's not as convenient as using logarithms nowadays. It comes from an algorithm by John von Neumann (1951) and uses only comparisons.

1. Let scale be 1/rate. Set highpart to 0.
2. Generate a uniform random variate in (0, scale), call it u.
3. Set val to u, and set accept to 1.
4. Generate a uniform random variate in (0, scale), call it v.
5. If u is greater than v, set u to v, then set accept to 1 minus accept, then go to step 4.
6. If accept is 1, return val + highpart.
7. Add scale to highpart and go to step 2.

REFERENCES:

• Von Neumann, J., "Various techniques used in connection with random digits", 1951.

---

If I understand your problem, and you can accept a finite number of PRNG's, you could follow an approach like:

• Create an array where every element is in your exponential distribution
• Generate a PRNG which is an integer index into the array. Return the element in the array at that index.

---

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.