29
$\begingroup$

You have one coin. You may flip it as many times as you want.

You want to generate a random number $r$ such that $a \leq r < b$ where $r,a,b\in \mathbb{Z}^+$.

Distribution of the numbers should be uniform.

It is easy if $b -a = 2^n$:

r = a + binary2dec(flip n times write 0 for heads and 1 for tails) 

What if $b-a \neq 2^n$?

$\endgroup$
2
  • $\begingroup$ Use Han-Hoshi algorithm - basically split the interval into two, use your random bit (coin flip) to randomly choose one of the two intervals, then repeat this process on the side you chose until you run out of bits. This will give you an interval uniformly distributed from a partition of the real line. The more flips you have, the more precise the interval. $\endgroup$
    – zenna
    Aug 16, 2015 at 16:38
  • $\begingroup$ I'd look at Arithmetic Coding. $\endgroup$
    – gnasher729
    Aug 2, 2022 at 14:48

10 Answers 10

14
$\begingroup$

What you're looking for is based on Rejection sampling or the accept-reject method (note that the Wiki page is a bit technical).

This method is useful in these kinds of situations: you want to pick some random object from a set (a random integer in the set $[a,b]$ in your case), but you don't know how to do that, but you can pick some random object from a larger set containing the first set (in your case, $[a, 2^k + a]$ for some $k$ such that $2^k + a \ge b$; this corresponds to $k$ coin flips).

In such a scenario, you just keep picking elements from the bigger set until you randomly pick an element in the smaller set. If your smaller set is big enough compared to your larger set (in your case, $[a, 2^k + a]$ contains at most twice as many integers as $[a,b]$, which is good enough), this is efficient.

An alternative example: suppose you want to pick a random point inside a circle with radius 1. Now, it isn't really easy to come up with a direct method for this. We turn to the accept-reject method: we sample points in a 1x1 square encompassing the circle, and test if the number we draw lies inside the circle.

$\endgroup$
2
  • 3
    $\begingroup$ Note that if we reject samples from $A$ in order to get a distribution on $B$, the expected number of iterations is $\frac{|A|}{|B|}$ (as we perform an experiment with geometric distribution). $\endgroup$
    – Raphael
    Mar 21, 2012 at 15:23
  • $\begingroup$ I recall seeing somewhere that this can't be done exactly unless the range is a power of 2 (as stands to reason, e.g. 1 / 3 has no terminating binary expansion). $\endgroup$
    – vonbrand
    Jan 25, 2013 at 14:22
7
$\begingroup$

(technicalities: the answer fits selection of number $a \le x < b$)

Since you are allowed to flip your coin as many times as you wish, you can get your probability as-close-as-you-wish to uniform by picking a fraction $r\in [0,1]$ (using binary radix: you flip the coin for each digit after the point) and multiply $r$ by $b-a$ to get a number between 0 and [b-a-1] (rounding down to an integer). Add this number to $a$ and you're done.

Example: say $b-a=3$. 1/3 in binary is 0.0101010101.... Then, if your flip is between 0 and 0.010101... your pick is $b$. if it is beween 0.010101.. and 0.10101010... your pick will be $a+1$, and otherwise it will be $a+2$.

If you flip your coin $t$ times then each number between $a$ and $b$ will be chosen with probability $\frac{1}{b-a}\pm 2^{-(t+1)}$.

$\endgroup$
6
  • 1
    $\begingroup$ This doesn't give a uniform distribution. For some applications (e.g. crypto, sometimes), this can be very bad. $\endgroup$ Mar 21, 2012 at 9:58
  • 3
    $\begingroup$ @Gilles: It can be fixed to give a perfectly uniform distribution by flipping until it's no longer possible for the result to change. This is the most efficient answer. $\endgroup$
    – Neil G
    Mar 21, 2012 at 11:24
  • $\begingroup$ @NeilG I know it can be fixed, but fixing it would be a crucial part of the answer. $\endgroup$ Mar 21, 2012 at 12:33
  • 2
    $\begingroup$ @Gilles: You're right. He could modify the answer to say that you can produce a perfectly uniform distribution if you flip while $\lfloor (b-a)(f + 2^{-t-1})\rfloor \ne \lfloor (b-a)(f - 2^{-t-1})\rfloor$. +1 from me for the best average case and worst case time. $\endgroup$
    – Neil G
    Mar 21, 2012 at 12:45
  • $\begingroup$ @NeilG, it can't be "fixed", as there is quite a sizeable set of integers which don't have a terminating binary fraction. $\endgroup$
    – vonbrand
    Mar 20, 2013 at 20:02
7
$\begingroup$

Choose a number in the next larger power of 2 range, and discard answers greater than $b-a$.

n = b-a;
N = round_to_next_larger_power_of_2(n)
while (1) {
  x = random(0 included to N excluded);
  if (x < n) break;
}
r = a + x;
$\endgroup$
3
  • 4
    $\begingroup$ And why does this work? $\endgroup$
    – Raphael
    Mar 21, 2012 at 13:20
  • $\begingroup$ @Raphael are you skeptical, or do you merely want the poster to explain in more detail ? $\endgroup$
    – Suresh
    Mar 21, 2012 at 15:30
  • 2
    $\begingroup$ @Suresh: The latter. The pseudo code could be polished a bit, but it implements what other answerers explain. Without justification, this answer is not worth much on its own. $\endgroup$
    – Raphael
    Mar 21, 2012 at 15:39
7
$\begingroup$

Nobody mentioned this, so let me formally prove that if $D$ is a domain whose size is not a power of two, then finitely many fair coin tosses aren't enough to generate a uniformly random member of $D$. Suppose you use $k$ coins to generate a member of $D$. For each $d \in D$, the probability that your algorithm generated $d$ is of the form $A/2^k$ for some integer $A$. The fundamental theorem of arithmetic shows that $A/2^k \neq 1/{|D|}$.

If you want to generate $n$ independent uniform samples of $D$, then the expected number of coin tosses you need (under the optimal algorithm) is roughly $n \log_2 |D|$. More generally, if you want to generate from a distribution of entropy $H$, you need roughly $nH$ random bits in expectation. An algorithm achieving this is arithmetic decoding, applied to an (ostensibly) infinite sequence of random bits.

$\endgroup$
3
  • $\begingroup$ However, it is possible to increase domain $D$ with dummy values until $\frac{A}{2^k}=\frac{1}{|D|}$, and then discard any dummy value outcomes - in which case, it is possible to generate a uniformly random member of the original domain with a finite number of coin tosses. $\endgroup$ Dec 13, 2019 at 1:54
  • $\begingroup$ That won’t be a bounded number of coin tosses. $\endgroup$ Dec 13, 2019 at 8:26
  • $\begingroup$ You're right - I see what you're saying. $\endgroup$ Dec 13, 2019 at 11:09
3
$\begingroup$

If b-a is not a power of 2 then you may have to flip many coins in order to get a result. You may even never get a result, but that is unlikely in the extreme.

Methods

The simplest method is to generate a number in [a, a+2^n), where 2^n >= b-a, until one happens to land in [a, b). You throw away a lot of entropy with this method.

A more expensive method allows you to keep all of the entropy, but becomes very expensive computationally as the number of coin flips / dice rolls increases. Intuitively it is like treating the coin flips as the digits of a binary number to the right of the decimal point, converting that number from base 2 to base a-b after, and returning the digits of that number as they become 'stuck'.

Example

The following code converts rolls of a fair n-sided die into rolls of a fair m-sided die (in your case n=2, m=a-b), with increasing marginal cost as the number of rolls increases. Note the need for a Rational number type with arbitrary precision. One nice property is that converting from n-sided to m-sided and back to n-sided will return the original stream, though perhaps delayed by a couple rolls due to digits having to become stuck.

public static IEnumerable<BigInteger> DigitConversion(this IEnumerable<BigInteger> inputStream, BigInteger modIn, BigInteger modOut) {
    //note: values are implicitly scaled so the first unfixed digit of the output ranges from 0 to 1
    Rational b = 0; //offset of the chosen range
    Rational d = 1; //size of the chosen range
    foreach (var r in inputStream) {
        //narrow the chosen range towards the real value represented by the input
        d /= modIn;
        b += d * r;
        //check for output digits that have become fixed
        while (true) {
            var i1 = (b * modOut).Floor();
            var i2 = ((b + d) * modOut).Floor(); //note: ideally b+d-epsilon, but another iteration makes that correction unnecessary
            if (i1 != i2) break; //digit became fixed?
            //fix the next output digit (rescale the range to make next digit range from 0 to 1)
            d *= modOut;
            b *= modOut;
            b -= i1;
            yield return i1;
        }
    }
}
$\endgroup$
1
  • $\begingroup$ "but that is unlikely in the extreme" -- the probability for this event is $0$; we say it "almost surely" does not happen. $\endgroup$
    – Raphael
    Mar 21, 2012 at 15:21
2
$\begingroup$

Generate a binary decimal. Instead of storing it explicitly, just keep track of the min and max possible values. Once those values lie within the same integer, return that integer. Sketch of the code is below.

(Edit) Fuller explanation: Say you want to generate a random integer from 1 to 3 inclusive with 1/3 probability each. We do this by generating a random binary decimal real, x in the range (0, 1). If x < 1/3, return 1, else if x < 2/3 return 2, else return 3. Instead of generating the digits for x explicitly, we just keep track of the minimum and maximum possible values of x. Initially, the min value of x is 0 and the max is 1. If you flip heads first then the first digit of x behind the decimal point (in binary) is 1. The minimum possible value of x (in binary) then becomes 0.100000 = 1/2 and the max is 0.111111111 = 1. Now if your next flip is tails x starts with 0.10. The min possible value is 0.1000000 = 1/2 and the max is 0.1011111 = 3/4. The min possible value of x is 1/2 so you know there's no chance of returning 1 since that requires x < 1/3. You can still return 2 if x ends up being 1/2 < x < 2/3 or 3 if 2/3 < x < 3/4. Now suppose the third flip is tails. Then x must start with 0.100. Min = 0.10000000 = 1/2 and max = 0.100111111 = 5/8. Now since 1/3 < 1/2 < 5/8 < 2/3 we know that x must fall in the interval (1/3, 2/3), so we can stop generating digits of x and just return 2.

The code does essentially this except instead of generating x between 0 and 1 it generates x between a and b, but the principle is the same.

def gen(a, b):
  min_possible = a
  max_possible = b

  while True:
    floor_min_possible = floor(min_possible)
    floor_max_possible = floor(max_possible)
    if max_possible.is_integer():
      floor_max_possible -= 1

    if floor_max_possible == floor_min_possible:
      return floor_max_possible

    mid = (min_possible + max_possible)/2
    if coin_flip():
      min_possible = mid
    else:
      max_possible = mid

Remark: I tested this code against the accept/reject method and both yield uniform distributions. This code requires less coin flips than accept reject except when b - a is close to the next power of 2. E.x. if you want to generate a = 0, b = 62 then accept/reject does better. I was able to prove that this code can use on average no more than 2 more coin flips than accept/reject. From my reading, it looks like Knuth and Yao (1976) gave a method to solve this problem and proved that their method is optimal in the expected number of coin flips. They further proved that expected number of flips must be greater than the Shannon entropy of the distribution. I could not find a copy of the text of the paper however and would be curious to see what their method is. (Update: just found an exposition of Knuth Yao 1976 here: http://www.nrbook.com/devroye/Devroye_files/chapter_fifteen_1.pdf but I have not yet read it). Someone also mentioned Han Hoshi in this thread which seems to be more general and solves it using a biased coin. See also http://paper.ijcsns.org/07_book/200909/20090930.pdf by Pae (2009) for a good discussion of the literature.

$\endgroup$
0
2
$\begingroup$

Theoretically optimal algorithm

Here is an improvement of the other answer I posted. The other answer has the advantage that it is easier to extend to the more general case of generating one discrete distribution from another. In fact, the other answer is a special case of the algorithm due to Han and Hoshi.

The algorithm I will describe here is based on Knuth and Yao (1976). In their paper, they also proved that this algorithm achieves the minimum possible expected number of coin flips.

To illustrate it, consider the Rejection sampling method described by other answers. As an example, suppose you want to generate one of 5 numbers uniformly [0, 4]. The next power of 2 is 8 so you flip the coin 3 times and generate a random number up to 8. If the number is 0 to 4 then you return it. Otherwise, you throw it out and generate another number up to 8 and try again until you succeed. But when you throw out the number, you just wasted some entropy. You can instead condition on the number you threw out to reduce the number of future coin flips you'll need in expectation. Concretely, once you generate the number [0, 7], if it is [0, 4], return. Otherwise, it's 5, 6, or 7, and you do something different in each case. If it's 5, flip the coin again and return either 0 or 1 based on the flip. If it's 6, flip the coin and return either 2 or 3. If it's 7, flip the coin; if it's heads, return 4, if it's tails start over.

The leftover entropy from our initial failed attempt gave us 3 cases (5, 6, or 7). If we just throw this out, we throw away log2(3) coin flips. We instead keep it, and combine it with the outcome of another flip to generate 6 possible cases (5H, 5T, 6H, 6T, 7H, 7T) which let's us immediately try again to generate a final answer with probability of success 5/6.

Here is the code:

# returns an int from [0, b)
def __gen(b):
  rand_num = 0
  num_choices = 1

  while True:
    num_choices *= 2
    rand_num *= 2
    if coin.flip():
      rand_num += 1

    if num_choices >= b:
      if rand_num < b:
        return rand_num
      num_choices -= b
      rand_num -= b

# returns an int from [a, b)
def gen(a, b):
  return a + __gen(b - a)
$\endgroup$
1
$\begingroup$

The simple answer?

If $b-a$ is not a power of 2, then after generating $r$, check whether it is in-range, and if not, generate it again.

The most likely time you will have to re-generate $r$ is when $b - a = 2^n + 1$ for some integer $n$, but even then it will have a greater than 50% chance of falling within range on each generation.

$\endgroup$
1
  • $\begingroup$ This doesn't seem complete. $\endgroup$ Mar 21, 2012 at 9:47
1
$\begingroup$

This is a proposed solution for the case when b - a does not equal 2^k. It is supposed to work in fixed number of steps (no need to throw away candidates which are outside your expected range).

However, I am not sure this correct. Please critique and help describe the exact non-uniformity in this random number generator (if any), and how to measure/quantify it.

Firstly convert to equivalent problem of generating uniformly distributed random numbers in range [0, z-1] where z = b - a.

Also, let m = 2^k be the smallest power of 2 >= z.

As per the solution above, we already have a uniformly distributed random number generator R(m) in range [0,m-1] (can be done by tossing k coins, one for each bit).

    Keep a random seed s and initialize with s = R(m).   

    function random [0, z-1] :
        x = R(m) + s 
        while x >= z:
            x -= z
        s = x
        return x

The while loop runs at most 3 times, giving next random number in fixed number of steps (best case = worst case).

See a test program for numbers [0,2] here: http://pastebin.com/zuDD2V6H

$\endgroup$
4
  • $\begingroup$ This is not uniform. Take $z=3$. You get $m=4$ and the probabilities are $1/2,1/4,1/4$. $\endgroup$ Jan 24, 2013 at 21:49
  • $\begingroup$ Please have a look at the pseudo code as well as the linked code more closely. It does emit 0, 1, and 2 almost with equal frequency... $\endgroup$
    – vpathak
    Jan 25, 2013 at 2:10
  • $\begingroup$ You're right, when you look at the outputs individually, they are uniform (the stationary probability is uniform). But they're not independent. If you have just output $0$, say, then the next output is zero with probability $1/2$ and one or two with probability $1/4$ each. $\endgroup$ Jan 25, 2013 at 4:08
  • $\begingroup$ You can replace the entire function with a single line: return s = (s + R(m)) % z; $\endgroup$ Jan 25, 2013 at 4:10
0
$\begingroup$

If $2^{n-1} < b - a ≤ 2^n$ then pick a number x, $0 ≤ x < 2^n$, using n coin throws to determine the highest, second highest etc. bit of x. If x < b-a then return a + x otherwise try again. As an optimisation, try again as soon as the throws so far tell you that $x < b - a$ is not possible. For example if b - a = 33, and the first coin is 1 = 32, and the second coin is 1 = 16, then x ≥ 48 so restart already.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.