views:

1300

answers:

6

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?

A: 

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.
GreenMatt
This approximate method works for difficult cases, but is not necessary this time out. An exponential distribution can be thrown exactly.
dmckee
A: 

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.

slebetman
Gets the wrong distribution.
dmckee
Yep. Originally I had **"Feel free to substitute with appropriate formula"** at the end which was accidentally deleted during an edit. The idea was to give an idea how to do this then google/read up what the appropriate formula would be\*. I'm not re-editing the answer because otherwise this comment would not make sense. Just read this comment before downvoting.
slebetman
Feel free to substitute this comment with one that is appropriate.
Charlie Salts
+2  A: 

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.

Ramashalanka
+7  A: 
Alok
Yep.. this is what I needed. After reading the wikipedia page more closely, the inversion method makes a lot of sense. For others reading your answer, there may be some confusion about the base of your log function. Technically, it should be base e, ie, ln( ).
Charlie Salts
Yes, my `log` is the natural log, although any base will do, just with a different lambda :-).
Alok
+3  A: 

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.

dmckee
Thorough answer, thank you.
Charlie Salts
A: 

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.

Grembo
Interesting idea. How do I control the λ value?
Charlie Salts
Oh - I misstated the process a bit. I actually created 1001 rv's on the interval (0,1000) and took 1000 differences. The result is exponential with mean 1.0 since the average difference is 1.0. To get another mean, just multiply the difference by the mean you want. BTW, I checked the results in @Risk to make sure the distribution was exponential with mean 1.0.
Grembo