Function returning random double with exponential distribution in range (a,b)

Question

I want to generate a random number from a to b. The problem is, the number has to be given with exponential distribution.

Here's my code:

public double getDouble(double low, double high)
        {
            double r;
            (..some stuff..)
            r = rand.NextDouble();
            if (r == 0) r += 0.00001;
            return (1 / -0.9) * Math.Log(1 - r) * (high - low) + low;
        }

The problem is that (1 / -0.9) * Math.Log(1 - r) is not between 0 and 1, so the result won't be between a and b. Can someone help? Thanks in advance!

Answer 1

I missunderstood your question in the first answer :) You are already using the inversion sampling.

To map a range into another range, there is a typical mathematical approach:

f(x) = (b-a)(x - min)/(max-min) + a

where

b = upper bound of target
a = lower bound of target
min = lower bound of source
max = upper bound of source
x = the value to map

(this is linear scaling, so the distribution would be preserved)

(You can verify: If you put in min for x, it results in a, if you put in max for x, you'll get b.)

The Problem now: The exponential distribution has a maximum value of inf. So, you cannot use this equation, because it always wold be whatever / inf + 0 - so 0. (Which makes sense mathematically, but ofc. does not fit your needs)

So, the ONLY correct answer is: There is no exponential distribution possible between two fixed numbers, cause you can't map [0,inf] -> [a,b]

Therefore you need some sort of trade-off, to make your result as exponential as possible.

I wrapped my head around different possibilities out of curiosity and I found that you simple can't beat maths on this :P

However, I did some test with Excel and 1.4 Million random records: I picked a random number as "limit" (10) and rounded the computed result to 1 decimal place. (0, 0.1, 0.2 and so on) This number I used to perform the linear transformation with an maximum of 10, ingoring any result greater than 1.

Out of 1.4 Million computations (generated it 10-20 times), only 7-10 random numbers greater than 1 have been generated:

(Probability density function, After mapping the values: Column 100 := 1, Column 0 := 0)

So:

• Map the values to [0,1], using the linear approach mentioned above, assume a maximum of 10 for the transformation.

• If you encounter a value > 1 after the transformation - just draw another random number, until the value is < 1.

• With only 7-10 occurences out of 1.4 Million tests, this should be close enough, since the re-drawn number will again be pseudo-exponential-distributed.

• If you want to build a spaceship, where navigation depends on perfectly exponential distributed numbers between 0 and 1 - don't do it, else you should be good.
• (If you want to cheat a bit: If you encounter a number > 1, just find the record that has the biggest variance (i.e. Max(occurrences < expected occurrences)) from it's expected value - then assume that value :P )

Answer 2

Since the support for the exponential distribution is 0 to infinity, regardless of the rate, I'm going to assume that you're asking for an exponential that's truncated below a and above b. Another way of expressing this would be an exponential random variable X conditioned on a <= X <= b.

You can derive the inversion algorithm for this by calculating the cumulative distribution function (CDF) of the truncated distribution as the integral from a to x of the density for your exponential. Scale the result by the area between a and b (which is F(b) - F(a) where F(x) is the CDF of the original exponential distribution) to make it a valid distribution with an area of 1. Set the derived CDF to U, a uniform(0,1) random number, and solve for X to get the inversion.

I don't program C#, but here's the result expressed in Ruby. It should translate pretty transparently.

def exp_in_range(a, b, rate = 1.0)
  exp_rate_a = Math.exp(-rate * a)
  return -Math.log(exp_rate_a - rand * (exp_rate_a - Math.exp(-rate * b))) / rate
end

I put a default rate of 1.0 since you didn't specify, but clearly you can override that. rand is Ruby's built-in uniform generator. I think the rest is pretty self-explanatory. I cranked out several test sets of 100k observations for a variety of (a,b) values, loaded the results into my favorite stats package, and the results are as expected.

Answer 3

The exponential distribution is not limited on the positive side, so values can go from 0 to inf. There are many ways to scale [0,infinity] to some finite interval, but the result would not be exponential distributed. If you just want a slice of the exponential distribution between a and b, you could simply draw r from [ra rb] such that -log(1-ra)=a and -log(1-rb)=b , i,e,

r=rand.NextDouble(); // assume this is between 0 and 1

ra=Math.Exp(-a)-1;
rb=Math.Exp(-b)-1;

rbound=ra+(rb-ra)*r;
return -Math.Log(1 - rbound);

Why check for r==0? I think you would want to check for the argument of the log to be >0, so check for r (or rbound int this case) ==1. Also not clear why the (1/-.9) factor??