Unbiased random number generator using a biased one

You have a biased random number generator that produces a 1 with a probability p and 0 with a probability (1-p). You do not know the value o开发者_如何转开发f p. Using this make an unbiased random number generator which produces 1 with a probability 0.5 and 0 with a probability 0.5.

Note: this problem is an exercise problem from Introduction to Algorithms by Cormen, Leiserson, Rivest, Stein.(clrs)

The events (p)(1-p) and (1-p)(p) are equiprobable. Taking them as 0 and 1 respectively and discarding the other two pairs of results you get an unbiased random generator.

In code this is done as easy as:

int UnbiasedRandom()
{
    int x, y;

    do
    {
        x = BiasedRandom();
        y = BiasedRandom();
    } while (x == y);

    return x;
}

The procedure to produce an unbiased coin from a biased one was first attributed to Von Neumann (a guy who has done enormous work in math and many related fields). The procedure is super simple:

• Toss the coin twice.
• If the results match, start over, forgetting both results.
• If the results differ, use the first result, forgetting the second.

The reason this algorithm works is because the probability of getting HT is p(1-p), which is the same as getting TH (1-p)p. Thus two events are equally likely.

I am also reading this book and it asks the expected running time. The probability that two tosses are not equal is z = 2*p*(1-p), therefore the expected running time is 1/z.

The previous example looks encouraging (after all, if you have a biased coin with a bias of p=0.99, you will need to throw your coin approximately 50 times, which is not that many). So you might think that this is an optimal algorithm. Sadly it is not.

Here is how it compares with the Shannon's theoretical bound (image is taken from this answer). It shows that the algorithm is good, but far from optimal.

You can come up with an improvement if you will consider that HHTT will be discarded by this algorithm, but in fact it has the same probability as TTHH. So you can also stop here and return H. The same is with HHHHTTTT and so on. Using these cases improves the expected running time, but are not making it theoretically optimal.

And in the end - python code:

import random

def biased(p):
    # create a biased coin
    return 1 if random.random() < p else 0

def unbiased_from_biased(p):
    n1, n2 = biased(p), biased(p)
    while n1 == n2:
        n1, n2 = biased(p), biased(p)

    return n1

p = random.random()
print p

tosses = [unbiased_from_biased(p) for i in xrange(1000)]
n_1 = sum(tosses)
n_2 = len(tosses) - n_1
print n_1, n_2

It is pretty self-explanatory, and here is an example result:

0.0973181652114
505 495

As you see, nonetheless we had a bias of 0.097, we got approximately the same number of 1 and 0

The trick attributed to von Neumann of getting two bits at a time, having 01 correspond to 0 and 10 to 1, and repeating for 00 or 11 has already come up. The expected value of bits you need to extract to get a single bit using this method is 1/p(1-p), which can get quite large if p is especially small or large, so it is worthwhile to ask whether the method can be improved, especially since it is evident that it throws away a lot of information (all 00 and 11 cases).

Googling for "von neumann trick biased" produced this paper that develops a better solution for the problem. The idea is that you still take bits two at a time, but if the first two attempts produce only 00s and 11s, you treat a pair of 0s as a single 0 and a pair of 1s as a single 1, and apply von Neumann's trick to these pairs. And if that doesn't work either, keep combining similarly at this level of pairs, and so on.

Further on, the paper develops this into generating multiple unbiased bits from the biased source, essentially using two different ways of generating bits from the bit-pairs, and giving a sketch that this is optimal in the sense that it produces exactly the number of bits that the original sequence had entropy in it.

You need to draw pairs of values from the RNG until you get a sequence of different values, i.e. zero followed by one or one followed by zero. You then take the first value (or last, doesn't matter) of that sequence. (i.e. Repeat as long as the pair drawn is either two zeros or two ones)

The math behind this is simple: a 0 then 1 sequence has the very same probability as a 1 then zero sequence. By always taking the first (or the last) element of this sequence as the output of your new RNG, we get an even chance to get a zero or a one.

Besides the von Neumann procedure given in other answers, there is a whole family of techniques, called randomness extraction (also known as debiasing, deskewing, or whitening), that serve to produce unbiased random bits from random numbers of unknown bias. They include Peres's (1992) iterated von Neumann procedure, as well as an "extractor tree" by Zhou and Bruck (2012). Both methods (and several others) are asymptotically optimal, that is, their efficiency (in terms of output bits per input) approaches the optimal limit as the number of inputs gets large (Pae 2018).

For example, the Peres extractor takes a list of bits (zeros and ones with the same bias) as input and is described as follows:

1. Create two empty lists named U and V. Then, while two or more bits remain in the input:
   • If the next two bits are 0/0, append 0 to U and 0 to V.
   • Otherwise, if those bits are 0/1, append 1 to U, then write a 0.
   • Otherwise, if those bits are 1/0, append 1 to U, then write a 1.
   • Otherwise, if those bits are 1/1, append 0 to U and 1 to V.
2. Run this algorithm recursively, reading from the bits placed in U.
3. Run this algorithm recursively, reading from the bits placed in V.

This is not to mention procedures that produce unbiased random bits from biased dice or other biased random numbers (not just biased bits); see, e.g., Camion (1974).

I discuss more on randomness extractors in a note on randomness extraction.

REFERENCES:

• Peres, Y., "Iterating von Neumann's procedure for extracting random bits", Annals of Statistics 1992,20,1, p. 590-597.
• Zhou, H. And Bruck, J., "Streaming algorithms for optimal generation of random bits", arXiv:1209.0730 [cs.IT], 2012.
• S. Pae, "Binarization Trees and Random Number Generation", arXiv:1602.06058v2 [cs.DS].
• Camion, Paul, "Unbiased die rolling with a biased die", North Carolina State University. Dept. Of Statistics, 1974.

Here's one way, probably not the most efficient. Chew through a bunch of random numbers until you get a sequence of the form [0..., 1, 0..., 1] (where 0... is one or more 0s). Count the number of 0s. If the first sequence is longer, generate a 0, if the second sequence is longer, generate a 1. (If they're the same, try again.)

This is like what HotBits does to generate random numbers from radioactive particle decay:

Since the time of any given decay is random, then the interval between two consecutive decays is also random. What we do, then, is measure a pair of these intervals, and emit a zero or one bit based on the relative length of the two intervals. If we measure the same interval for the two decays, we discard the measurement and try again

HotBits: How It Works

I'm just explaining the already proposed solutions with some running proof. This solution will be unbiased, no matter how many times we change the probability. In a head n tail toss, the exclusivity of consecutive head n tail or tail n head is always unbiased.

import random

def biased_toss(probability):
    if random.random() > probability:
        return 1
    else:
        return 0
    
def unbiased_toss(probability):
    x = biased_toss(probability)
    y = biased_toss(probability)
    while x == y:
        x = biased_toss(probability)
        y = biased_toss(probability)
    else:
        return x

# results with contain counts of heads '0' and tails '1'
results = {'0':0, '1':0}
for i in range(1000):
    # on every call we are changing the probability
    p = random.random()
    results[str(unbiased_toss(p))] += 1

# it still return unbiased result
print(results)