Generating evenly distributed bits, using approximation

Question

I'm trying to generate 0 or 1 with 50/50 chance of any using random.uniform instead of random.getrandbits.

Here's what I have

0 if random.uniform(0, 1e-323) == 0.0 else 1

But if I run this long enough, the average is ~70% to generate 1. As seem here:

sum(0 if random.uniform(0, 1e-323) == 0.0 
    else 1 
    for _ in xrange(1000)) / 1000.0  # --> 0.737

If I change it to 1e-324 it will always be 0. And if I change it to 1e-322, the average will be ~%90.

I made a dirty program that will try to find the sweet spot between 1e-322 and 1e-324, by dividing and multiplying it several times:

v = 1e-323
n_runs = 100000
target = n_runs/2

result = 0
while True:
    result = sum(0 if random.uniform(0, v) == 0.0 else 1 for _ in xrange(n_runs))

    if result > target:
        v /= 1.5
    elif result < target:
        v *= 1.5 / 1.4
    else:
        break

print v

This end ups with 4.94065645841e-324

But it still will be wrong if I ran it enough times.

Is there I way to find this number without the dirty script I wrote? I know that Python has a intern min float value, show in sys.float_info.min, which in my PC is 2.22507385851e-308. But I don't see how to use it to solve this problem.

Sorry if this feels more like a puzzle than a proper question, but I'm not able to answer it myself.

Answer 1

I know that Python has a intern min float value, show in sys.float_info.min, which in my PC is 2.22507385851e-308. But I don't see how to use it to solve this problem.

2.22507385851e-308 is not the smallest positive float value, it is the smallest positive normalized float value. The smallest positive float value is 2^-52 times that, that is, near 5e-324.

2^-52 is called the “machine epsilon” and it is usual to call the “min” of a floating-point type a value that is nether that which is least of all comparable values (that is -inf), nor the least of finite values (that is -max), nor the least of positive values.

Then, the next problem you face is that random.uniform is not uniform to that level. It probably works ok when you pass it a normalized number, but if you pass it the smallest positive representable float number, the computation it does with it internally may be very approximative and lead it to behave differently than the documentation says. Although it appears to work surprisingly ok according to the results of your “dirty script”.

Answer 2

Here's the random.uniform implementation, according to the source:

from os import urandom as _urandom

BPF = 53        # Number of bits in a float
RECIP_BPF = 2**-BPF

def uniform(self, a, b):
    "Get a random number in the range [a, b) or [a, b] depending on rounding."
     return a + (b-a) * self.random()

def random(self):
     """Get the next random number in the range [0.0, 1.0)."""
     return (int.from_bytes(_urandom(7), 'big') >> 3) * RECIP_BPF

So, your problem boils down to finding a number b that will give 0 when multiplied by a number less than 0.5 and another result when multiplied by a number larger than 0.5. I've found out that, on my machine, that number is 5e-324.

To test it, I've made the following script:

from random import uniform

def test():
    runs = 1000000
    results = [0, 0]
    for i in range(runs):
        if uniform(0, 5e-324) == 0:
            results[0] += 1
        else:
            results[1] += 1
    print(results)

Which returned results consistent with a 50% probability:

>>> test()
[499982, 500018]
>>> test()
[499528, 500472]
>>> test()
[500307, 499693]