Exponentially distributed random generator (log fu

Question

I really need help as I am stuck at the begining of the code.

I am asked to create a function to investigate the exponential distribution on histogram. The function is x = −log(1−y)/λ. λ is a constant and I referred to that as lamdr in the code and simply gave it 10. I gave N (the number of random numbers) 10 and ran the code yet the results and the generated random numbers gave me totally different results; below you can find the code, I don't know what went wrong, hope you guys can help me!! (I use python 2)

import random
import math

N = raw_input('How many random numbers you request?: ')
N = int(N)
lamdr = raw_input('Enter a value:')
lamdr = int(lamdr)

def exprand(lamdr):
    y = []
    for i in range(N):
        y.append(random.uniform(0,1))
    return y
y = exprand(lamdr)
print 'Randomly generated numbers:', (y)

x = []
for w in y:
    x.append((math.log((1 - w) / lamdr)) * -1)
print 'Results:', x

Answer 1

WHat @pjs wrote is true to a point. While statement mathematically, if y has a uniform(0,1) distribution, so does 1-y appears to be correct, proposal to replace code with -math.log(random.random()) / lambdr is just wrong. Why? Because Python random module provide U(0,1) in the range [0,1) (as mentioned here), thus making such replacement non-equivalent.

In more layman term, if your U(0,1) is actually generating numbers in the [0,1) range, then code

import random
def exprand(lambda):
    return -math.log(1.0 - random.random()) / lambda

is correct, but code

import random
def exprand(lambda):
    return -math.log(random.random()) / lambda

is wrong, it will sometimes generate NaN/exception, as log(0) will be called

Answer 2

After viewing the code you provided, it looks like you have the pieces you need but you're not putting them together.

You were asked to write function exprand(lambdr) using the specified formula. Python already provides a function called random.expovariate(lambd) for generating exponentials, but what the heck, we can still make our own. Your formula requires a "random" value for y which has a uniform distribution between zero and one. The documentation for the random module tells us that random.random() will give us a uniform(0,1) distribution. So all we have to do is replace y in the formula with that function call, and we're in business:

def exprand(lambdr):
    return -math.log(1.0 - random.random()) / lambdr

An historical note: Mathematically, if y has a uniform(0,1) distribution, then so does 1-y. Implementations of the algorithm dating back to the 1950's would often leverage this fact to simplify the calculation to -math.log(random.random()) / lambdr. Mathematically this gives distributionally correct results since P{X = c} = 0 for any continuous random variable X and constant c, but computationally it will blow up in Python for the 1 in 2⁶⁴ occurrence where you get a zero from random.random(). One historical basis for doing this was that when computers were many orders of magnitude slower than now, ditching the one additional arithmetic operation was considered worth the minuscule risk. Another was that Prime Modulus Multiplicative PRNGs, which were popular at the time, never yield a zero. These days it's primarily of historical interest, and an interesting example of where math and computing sometimes diverge.

Back to the problem at hand. Now you just have to call that function N times and store the results somewhere. Likely candidates to do so are loops or list comprehensions. Here's an example of the latter:

abuncha_exponentials = [exprand(0.2) for _ in range(5)]

That will create a list of 5 exponentials with λ=0.2. Replace 0.2 and 5 with suitable values provided by the user, and you're in business. Print the list, make a histogram, use it as input to something else...

Replacing exporand with expovariate in the list comprehension should produce equivalent results using Python's built-in exponential generator. That's the beauty of functions as an abstraction, once somebody writes them you can just use them to your heart's content.

Note that because of the use of randomness, this will give different results every time you run it unless you "seed" the random generator to the same value each time.