Why do people say there is modulo bias when using

Question

I have seen this question asked a lot but never seen a true concrete answer to it. So I am going to post one here which will hopefully help people understand why exactly there is "modulo bias" when using a random number generator, like rand() in C++.

Answer 1

I just wrote a code for Von Neumann's Unbiased Coin Flip Method, that should theoretically eliminate any bias in the random number generation process. More info can be found at (http://en.wikipedia.org/wiki/Fair_coin)

int unbiased_random_bit() {    
    int x1, x2, prev;
    prev = 2;
    x1 = rand() % 2;
    x2 = rand() % 2;

    for (;; x1 = rand() % 2, x2 = rand() % 2)
    {
        if (x1 ^ x2)      // 01 -> 1, or 10 -> 0.
        {
            return x2;        
        }
        else if (x1 & x2)
        {
            if (!prev)    // 0011
                return 1;
            else
                prev = 1; // 1111 -> continue, bias unresolved
        }
        else
        {
            if (prev == 1)// 1100
                return 0;
            else          // 0000 -> continue, bias unresolved
                prev = 0;
        }
    }
}

Answer 2

As the accepted answer indicates, "modulo bias" has its roots in the low value of RAND_MAX. He uses an extremely small value of RAND_MAX (10) to show that if RAND_MAX were 10, then you tried to generate a number between 0 and 2 using %, the following outcomes would result:

rand() % 3   // if RAND_MAX were only 10, gives
output of rand()   |   rand()%3
0                  |   0
1                  |   1
2                  |   2
3                  |   0
4                  |   1
5                  |   2
6                  |   0
7                  |   1
8                  |   2
9                  |   0

So there are 4 outputs of 0's (4/10 chance) and only 3 outputs of 1 and 2 (3/10 chances each).

So it's biased. The lower numbers have a better chance of coming out.

But that only shows up so obviously when RAND_MAX is small. Or more specifically, when the number your are modding by is large compared to RAND_MAX.

A much better solution than looping (which is insanely inefficient and shouldn't even be suggested) is to use a PRNG with a much larger output range. The Mersenne Twister algorithm has a maximum output of 4,294,967,295. As such doing MersenneTwister::genrand_int32() % 10 for all intents and purposes, will be equally distributed and the modulo bias effect will all but disappear.

Answer 3

Keep selecting a random is a good way to remove the bias.

Update

We could make the code fast if we search for an x in range divisible by n.

// Assumptions
// rand() in [0, RAND_MAX]
// n in (0, RAND_MAX]

int x; 

// Keep searching for an x in a range divisible by n 
do {
    x = rand();
} while (x >= RAND_MAX - (RAND_MAX % n)) 

x %= n;

The above loop should be very fast, say 1 iteration on average.

Answer 4

With a RAND_MAX value of 3 (in reality it should be much higher than that but the bias would still exist) it makes sense from these calculations that there is a bias:

1 % 2 = 1 2 % 2 = 0 3 % 2 = 1 random_between(1, 3) % 2 = more likely a 1

In this case, the % 2 is what you shouldn't do when you want a random number between 0 and 1. You could get a random number between 0 and 2 by doing % 3 though, because in this case: RAND_MAX is a multiple of 3.

Another method

There is much simpler but to add to other answers, here is my solution to get a random number between 0 and n - 1, so n different possibilities, without bias.

• the number of bits (not bytes) needed to encode the number of possibilities is the number of bits of random data you'll need
• encode the number from random bits
• if this number is >= n, restart (no modulo).

Really random data is not easy to obtain, so why use more bits than needed.

Below is an example in Smalltalk, using a cache of bits from a pseudo-random number generator. I'm no security expert so use at your own risk.

next: n

    | bitSize r from to |
    n < 0 ifTrue: [^0 - (self next: 0 - n)].
    n = 0 ifTrue: [^nil].
    n = 1 ifTrue: [^0].
    cache isNil ifTrue: [cache := OrderedCollection new].
    cache size < (self randmax highBit) ifTrue: [
        Security.DSSRandom default next asByteArray do: [ :byte |
            (1 to: 8) do: [ :i |    cache add: (byte bitAt: i)]
        ]
    ].
    r := 0.
    bitSize := n highBit.
    to := cache size.
    from := to - bitSize + 1.
    (from to: to) do: [ :i |
        r := r bitAt: i - from + 1 put: (cache at: i)
    ].
    cache removeFrom: from to: to.
    r >= n ifTrue: [^self next: n].
    ^r

Answer 5

Definition

Modulo Bias is the inherent bias in using modulo arithmetic to reduce an output set to a subset of the input set. In general, a bias exists whenever the mapping between the input and output set is not equally distributed, as in the case of using modulo arithmetic when the size of the output set is not a divisor of the size of the input set.

This bias is particularly hard to avoid in computing, where numbers are represented as strings of bits: 0s and 1s. Finding truly random sources of randomness is also extremely difficult, but is beyond the scope of this discussion. For the remainder of this answer, assume that there exists an unlimited source of truly random bits.

Problem Example

Let's consider simulating a die roll (0 to 5) using these random bits. There are 6 possibilities, so we need enough bits to represent the number 6, which is 3 bits. Unfortunately, 3 random bits yields 8 possible outcomes:

000 = 0, 001 = 1, 010 = 2, 011 = 3
100 = 4, 101 = 5, 110 = 6, 111 = 7

We can reduce the size of the outcome set to exactly 6 by taking the value modulo 6, however this presents the modulo bias problem: 110 yields a 0, and 111 yields a 1. This die is loaded.

Potential Solutions

Approach 0:

Rather than rely on random bits, in theory one could hire a small army to roll dice all day and record the results in a database, and then use each result only once. This is about as practical as it sounds, and more than likely would not yield truly random results anyway (pun intended).

Approach 1:

Instead of using the modulus, a naive but mathematically correct solution is to discard results that yield 110 and 111 and simply try again with 3 new bits. Unfortunately, this means that there is a 25% chance on each roll that a re-roll will be required, including each of the re-rolls themselves. This is clearly impractical for all but the most trivial of uses.

Approach 2:

Use more bits: instead of 3 bits, use 4. This yield 16 possible outcomes. Of course, re-rolling anytime the result is greater than 5 makes things worse (10/16 = 62.5%) so that alone won't help.

Notice that 2 * 6 = 12 < 16, so we can safely take any outcome less than 12 and reduce that modulo 6 to evenly distribute the outcomes. The other 4 outcomes must be discarded, and then re-rolled as in the previous approach.

Sounds good at first, but let's check the math:

4 discarded results / 16 possibilities = 25%

In this case, 1 extra bit didn't help at all!

That result is unfortunate, but let's try again with 5 bits:

32 % 6 = 2 discarded results; and
2 discarded results / 32 possibilities = 6.25%

A definite improvement, but not good enough in many practical cases. The good news is, adding more bits will never increase the chances of needing to discard and re-roll. This holds not just for dice, but in all cases.

As demonstrated however, adding an 1 extra bit may not change anything. In fact if we increase our roll to 6 bits, the probability remains 6.25%.

This begs 2 additional questions:

1. If we add enough bits, is there a guarantee that the probability of a discard will diminish?
2. How many bits are enough in the general case?

General Solution

Thankfully the answer to the first question is yes. The problem with 6 is that 2^x mod 6 flips between 2 and 4 which coincidentally are a multiple of 2 from each other, so that for an even x > 1,

[2^x mod 6] / 2^x == [2^(x+1) mod 6] / 2^(x+1)

Thus 6 is an exception rather than the rule. It is possible to find larger moduli that yield consecutive powers of 2 in the same way, but eventually this must wrap around, and the probability of a discard will be reduced.

Without offering further proof, in general using double the number of bits required will provide a smaller, usually insignificant, chance of a discard.

Proof of Concept

Here is an example program that uses OpenSSL's libcrypo to supply random bytes. When compiling, be sure to link to the library with -lcrypto which most everyone should have available.

#include <iostream>
#include <assert.h>
#include <limits>
#include <openssl/rand.h>

volatile uint32_t dummy;
uint64_t discardCount;

uint32_t uniformRandomUint32(uint32_t upperBound)
{
    assert(RAND_status() == 1);
    uint64_t discard = (std::numeric_limits<uint64_t>::max() - upperBound) % upperBound;
    uint64_t randomPool = RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));

    while(randomPool > (std::numeric_limits<uint64_t>::max() - discard)) {
        RAND_bytes((uint8_t*)(&randomPool), sizeof(randomPool));
        ++discardCount;
    }

    return randomPool % upperBound;
}

int main() {
    discardCount = 0;

    const uint32_t MODULUS = (1ul << 31)-1;
    const uint32_t ROLLS = 10000000;

    for(uint32_t i = 0; i < ROLLS; ++i) {
        dummy = uniformRandomUint32(MODULUS);
    }
    std::cout << "Discard count = " << discardCount << std::endl;
}

I encourage playing with the MODULUS and ROLLS values to see how many re-rolls actually happen under most conditions. A sceptical person may also wish to save the computed values to file and verify the distribution appears normal.

Answer 6

Mark's Solution (The accepted solution) is Nearly Perfect.

int x;

do {
    x = rand();
} while (x >= (RAND_MAX - RAND_MAX % n));

x %= n;

edited Mar 25 '16 at 23:16

Mark Amery 39k21170211

However, it has a caveat which discards 1 valid set of outcomes in any scenario where RAND_MAX (RM) is 1 less than a multiple of N (Where N = the Number of possible valid outcomes).

ie, When the 'count of values discarded' (D) is equal to N, then they are actually a valid set (V), not an invalid set (I).

Using Mark's Solution, Values are Discarded when: X => RM - RM % N

EG: 

Ran Max Value (RM) = 255
Valid Outcome (N) = 4

When X => 252, Discarded values for X are: 252, 253, 254, 255

So, if Random Value Selected (X) = {252, 253, 254, 255}

Number of discarded Values (I) = RM % N + 1 == N

 IE:

 I = RM % N + 1
 I = 255 % 4 + 1
 I = 3 + 1
 I = 4

   X => ( RM - RM % N )
 255 => (255 - 255 % 4) 
 255 => (255 - 3)
 255 => (252)

 Discard Returns $True

As you can see in the example above, when the value of X (the random number we get from the initial function) is 252, 253, 254, or 255 we would discard it even though these four values comprise a valid set of returned values.

IE: When the count of the values Discarded (I) = N (The number of valid outcomes) then a Valid set of return values will be discarded by the original function.

If we describe the difference between the values N and RM as D, ie:

D = (RM - N)

Then as the value of D becomes smaller, the Percentage of unneeded re-rolls due to this method increases at each natural multiplicative. (When RAND_MAX is NOT equal to a Prime Number this is of valid concern)

EG:

RM=255 , N=2 Then: D = 253, Lost percentage = 0.78125%

RM=255 , N=4 Then: D = 251, Lost percentage = 1.5625%
RM=255 , N=8 Then: D = 247, Lost percentage = 3.125%
RM=255 , N=16 Then: D = 239, Lost percentage = 6.25%
RM=255 , N=32 Then: D = 223, Lost percentage = 12.5%
RM=255 , N=64 Then: D = 191, Lost percentage = 25%
RM=255 , N= 128 Then D = 127, Lost percentage = 50%

Since the percentage of Rerolls needed increases the closer N comes to RM, this can be of valid concern at many different values depending on the constraints of the system running he code and the values being looked for.

To negate this we can make a simple amendment As shown here:

 int x;

 do {
     x = rand();
 } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

 x %= n;

This provides a more general version of the formula which accounts for the additional peculiarities of using modulus to define your max values.

Examples of using a small value for RAND_MAX which is a multiplicative of N.

Mark'original Version:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X >= (RAND_MAX - ( RAND_MAX % n ) )
When X >= 2 the value will be discarded, even though the set is valid.

Generalized Version 1:

RAND_MAX = 3, n = 2, Values in RAND_MAX = 0,1,2,3, Valid Sets = 0,1 and 2,3.
When X > (RAND_MAX - ( ( RAND_MAX % n  ) + 1 ) % n )
When X > 3 the value would be discarded, but this is not a vlue in the set RAND_MAX so there will be no discard.

Additionally, in the case where N should be the number of values in RAND_MAX; in this case, you could set N = RAND_MAX +1, unless RAND_MAX = INT_MAX.

Loop-wise you could just use N = 1, and any value of X will be accepted, however, and put an IF statement in for your final multiplier. But perhaps you have code that may have a valid reason to return a 1 when the function is called with n = 1...

So it may be better to use 0, which would normally provide a Div 0 Error, when you wish to have n = RAND_MAX+1

Generalized Version 2:

int x;

if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

    x %= n;
} else {
    x = rand();
}

Both of these solutions resolve the issue with needlessly discarded valid results which will occur when RM+1 is a product of n.

The second version also covers the edge case scenario when you need n to equal the total possible set of values contained in RAND_MAX.

The modified approach in both is the same and allows for a more general solution to the need of providing valid random numbers and minimizing discarded values.

To reiterate:

The Basic General Solution which extends mark's example:

 int x;

 do {
     x = rand();
 } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

 x %= n;

The Extended General Solution which Allows one additional scenario of RAND_MAX+1 = n:

int x;

if n != 0 {
    do {
        x = rand();
    } while (x > (RAND_MAX - ( ( ( RAND_MAX % n ) + 1 ) % n) );

    x %= n;
} else {
    x = rand();
}