Probability and Random Numbers

Date: 02/28/2004 at 10:23:49
From: GVK
Subject: Probability and Random Numbers

Let me start with an example first.  If I have a six-faced die, I can 
be sure that rolling it will produce 6 events of equal probability 
(1/6 each).  Suppose I blacken one of the faces and roll the die.  I 
will ignore the result and roll the die again if the blackened face 
shows up.  Can I be sure that the other five will show up with a 
probability of 1/5 each? 

The original problem asks if I can use a coin (which has 2 events of 
equal probability) to devise three events with 1/3rd probability each?  
I find devising 4 events easy; toss the coin twice and interpret the 
results 00, 01, 10, 11 (ie HH, HT, TH, TT) as the 4 events.  But, is 
it okay to disregard any one result (say 00), similar to blackening 
the die above, and claim that the other three events are of equal 
probability?  I find it hard to believe that a coin can be used to 
have three events of equal probability.

I tried devising various versions of the problem, like using the 
function rand5() to generate a function rand7(), and the blackened 
six-faced die, etc.

Date: 02/29/2004 at 09:15:20
From: Doctor Jacques
Subject: Re: Probability and Random Numbers

Hi GVK,

The method does work.  However, there is a catch; the number of coin 
tosses is theoretically unbounded.  Although it is quite unlikely, you 
can have a very long run of "00" results.  This means that you cannot 
guarantee that you will generate the next number in a predetermined 
amount of time.  There is no method that would allow you to generate 
those events with a finite number of bits, since 3 never divides any 
power of 2 (another way to see this is to notice that the binary 
expansion of 1/3 is a repeating infinite string of bits).

As the probability of failure is 1/4, the expected number of coin 
tosses is:

  2*(1 + 1/4 + (1/4)^2 ...) = 8/3 = 2.66...

This is rather wasteful, as the theoretical number of bits required 
to generate an event with probability 1/3 is:

  log(3) / log(2) = 1.585

If you want to generate a single event, you cannot do better.  If you 
want to produce a large number of such events, with an exact 
probability of 1/3, and without wasting too many bits, you can still 
do better.

The first question would be "why bother?".  There are several reasons
for doing this, beyond the fact that this is fun!

In some applications, such as those related to high-security
cryptography, it is really important to have exact probabilities and
truly random numbers.  The numbers returned by the random() routines
found in programming languages are not random at all; they are
produced by a deterministic algorithm.  They are called random because
their distributions are statistically similar to true random numbers.

When really random numbers are required, we use external devices
based, for example, on radioactive decay or thermal noise--these 
processes cannot be reproduced.  If you are interested in this, you
could read:

     http://www.fourmilab.ch/hotbits/ 

However, those devices usually have a very low rate; if you want to 
produce large amounts of numbers in a short time, you may want to use 
as few random bits as possible.

Let us look at a method for generating random integers between 0 and 
(n - 1), with an exact uniform distribution.  This is similar to the 
random(n) function.  However, that function, as it is implemented, is 
not exact; it generates integers in a large range (for example,
between 0 and 65535), and returns the integer part of (n*x)/65536.  If 
n is not a power of 2, some numbers will have a slightly larger 
probability than others, although the difference will be small if n 
is small.  In this case, we will generate numbers between 0 and 
(n - 1) with an exact probability of 1/n.

We assume that we have available a function NextBit(), that returns a 
single random bit (and we assume that calls to that function are 
expensive).  We assume also that we have to generate many numbers, so 
that only the average cost matters.

A first idea would be to collect many random bits, and to generate 
random numbers in batches.  For example, for a probability of 1/3, we 
could collect 8 bits at a time to produce a number between 0 and 255. 
As 3^5 = 243 < 256, we could discard the result if it is greater than 
242; otherwise, we will produce five numbers between 0 and 2. The 
probability of failure is 13/256, so the average number of 
experiments is:

  1 + 13/256 + (13/256)^2 ... = 256/243 = 1.053...

In each experiment, we consume 8 bits and produce 5 numbers, so the 
average number of bits per number is:

  1.053 * 8/5 = 1.686

which is already much closer to the theoretical figure of 1.585.

We can still do better.  The idea is that, in case of failure, the 
difference (x - 243) is a random number between 0 and 12, which is 
equivalent to log(13)/log(2) = 3.7 bits, and we could try to take 
advantage of it.

The following algorithm uses a similar (but slightly different) idea. 
It uses two internal integer variables, m and r, which are not reset 
at the beginning of the algorithm (in C, you would declare them 
as "static").  Initially, m = 1 and r = 0.  We also have a parameter
N, which is a large integer such that 2N can still be represented 
exactly in the computer.  As said before, n is the modulus of the 
numbers you want to produce (they will be between 0 and (n - 1)), we 
must have n < N, and we have a function NextBit() that returns a 
truly random bit.

1. WHILE m < N DO r : = 2*r + NextBit(); m = 2*m; (r is a random
   variable of modulus m)

2. Divide m by n : m = n*q + b

3. IF r >= n*q, let m : = m - n*q, r : = r - n*q (r is still a random
   variable of modulus m), and go to step 1.

4. Otherwise, let x : = r mod n, r : = [r/n], and m : = q, and return 
   x.

In step 4, as r is between 0 and (n*q - 1), we can consider r as a 
random variable of modulus n*q.  As this is divisble by n, (r mod n) 
will be uniformly distributed, and the quotient [r/n] will be 
uniformly distributed between 0 and n - 1.

Note that the running time is still unbounded, since you can return 
to step 1 arbitrarily often, but, in practice, this will give you an 
exact distribution while consuming a number of bits quite close to 
the theoretical minimum.

For example, assume that we want to generate decimal digits (n = 10) 
and we let N = 256.  A typical run would be as follows:

* We start with m = 1, r = 0, and get 8 bits. Assume they are
  1, 1, 1, 1, 1, 1, 0, 1.
  At the end of step 1, we have m = 256, r = 253.

* We divide 256 by 10, giving q = 25, n*q = 250.

* As 253 >= 250, we are out of luck.  The excess, r - n*q = 253 - 250
  = 3, is a random variable of modulus 256 - 250 = 6. We let
  therefore m = 6 and r = 3 (note that we have saved more than two
  bits, instead of throwing them away).

* We start again at step 1.  This time, we need to generate 6 bits, in
  order to have m >= 256.  Assume the process runs as follows
  (starting with m = 6, r = 3:

   NextBit() |  m  |  r
  -----------+-----+----
      0      |  12 |   6
      1      |  24 |  13
      1      |  48 |  27
      0      |  96 |  54
      0      | 192 | 108
      1      | 384 | 217

* Step 2 gives q = 38, n*q = 380

* As 217 < 380, we will return (217 mod 10) = 7.  As we were in the
  21st decade out of 38, we let r = 21, m = 38 to be used for the
  next call.  This means that we already have more than 5 bits in
  store.

We have used 14 bits to generate the first number, but this is merely 
the "startup cost".  As m = 38, the next call will need only 3 new 
bits in step 1.  In the long run, the average number of bits per 
number will be close to the theoretical value log(10)/log(2) = 3.32. 
The probability of failure at step 3 is at most 9/m <= 9/256.

A nice feature of this algorithm is that n is used only inside the 
routine.  You can use a different n for each call, without re-
initializing the variables m and r.

Does this help?  Write back if you'd like to talk about this some 
more, or if you have any other questions.

- Doctor Jacques, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 03/10/2004 at 10:43:15
From: GVK
Subject: Probability and Random Numbers

Hello Dr. Jacques,

Thanks a lot for your reply.  I have just one question.  In your 
reply, you wrote:

   since 3 never divides any power of 2 (another way to see this is to
   notice that the binary expansion of 1/3 is a repeating infinite
   string of bits).

Are the two equivalent?  5 never divides any power of 2 either, 
however 1/5 is not an infinite sequence of bits.  Is it possible to 
get 5 events of 1/5 probability using a coin?

Date: 03/12/2004 at 02:04:45
From: Doctor Jacques
Subject: Re: Probability and Random Numbers

Hi again GVK,

The two are equivalent.  In base N, the N-ary expansion of 1/k is 
infinite (and repeating) unless all prime factors of k are also prime 
factors of N, which amounts to saying that k divides some power of N. 
If N = 2, the only allowed prime factor is 2, so any fraction whose 
denominator (in lowest terms) is not a power of 2 will have an 
infinite binary expansion.

For your specific example, the binary expansion of 1/5 is:

  1/5 = 0. 0011 0011.... = (3/16) + (3/16)^2 + ...

and this is a repeating binary number.

It is indeed quite possible to generate events with probability 1/5 
by tossing coins--you can use the algorithm I gave you with n = 5 
and N >= 8).

By the way, note that that algorithm also works in the opposite 
direction; you can also use it to convert "base n" events to "base 
m" events for any different m and n.  For example, you can use it to 
generate random binary strings using six-faced dice; you must 
essentially modify step 1 and the choice of N.

- Doctor Jacques, The Math Forum
  http://mathforum.org/dr.math/