15 views (last 30 days)

Show older comments

The Birthday Paradox or problem asks for the probability that in a room of n people, 2 or more have the same birthday (not date), assuming all years have N = 365 days. It is called a paradox because most people are surprised by the answer when there are (say) 30 people in the room.

Many applications require long sequences (or large vectors) of random numbers. In MATLAB these are supplied by the simple statement r = rand(n,1). This statement fills the vector r(n,1) with double precision random numbers uniformly distributed on (0,1).

Here is my question:

What is the probability that duplicates occur in r = rand(n,1), as n gets large?

I have an answer to this question but I would like to see what others come up with so that I can check the validity of my answer.

Derek O'Connor
on 5 Feb 2011

Walter Roberson
on 6 Feb 2011

Derek:

>> num2hex(2^-53)

ans =

3ca0000000000000

>> num2hex(1-2^(-53))

ans =

3fefffffffffffff

That is 53 different exponents that are used over the range of random numbers you have documented.

The only integral multiple of 2^(-53) that has an exponent of 3ca is 2^(-53) itself, with an all-0 mantissa. 2*2^(-53) has an exponent of 3cb and an all-0 mantissa. 3*2^(-53) has an exponent of 3cb and a mantissa whose first bit is set. 4*2^(-53) has an exponent of 3cc and a mantissa of 0.

If we follow this further than it becomes clear from the pigeon-hole principle that not all of the values starting with 3cf can be represented.

Example:

>> num2hex(hex2num('3fe0ffffffffffff') / 2^(-53))

ans =

4330ffffffffffff

It follows that 3fe0ffffffffffff is not an exact integral multiple of 2^(-53) and thus will not be randomly generated. (it is the representable number just below 0.53125 exactly)

Any argument made on the basis of equi-probable bit patterns after a constant exponent are suspect. One _can_, however, make arguments about taking integer multiples of 2^(-53) and then normalizing the numbers.

IEEE 754 does not actually store 53 bits of mantissa: it uses 52 bits of mantissa plus the "hidden bit" implied by the exponent.

Richard Willey
on 7 Feb 2011

Hi Derek

Your question about the birthday problem is best viewed as a special case of a more general topic: How should measure the statistical properties of a random number generator?

Rather than using this formulation of the birthday problem to reinvent the wheel, I recommend looking at established literature on this topic. A lot of very qualified experts have spent a lot of time and effort developing criteria to measure the quality of random number generators. The “Diehard” test suite is one such test (Diehard uses a reformulated version of the birthday paradox as one part of the test suite). Alternatively, there is another test suite that uses a combination of

- A monobit test (testing the distribution of ones and zeros in a sequence)
- A “Poker” test (a modified chi-square test)
- A Longruns test (checking whether there are runs of length 34 or greater within a 20,000 bit sequence)
- An autocorrelation test

Most of MATLAB’s code can be inspected by users. If you have questions about the specific implementation of a given algorithm – for example, what version of Twister does MathWorks use – you typically have the option to inspect the function and see how the code was written.

Regards,

Richard

Richard Willey
on 9 Feb 2011

Hi Derek

Sorry that I misconstrued your original question. Here’s one last quick comment:

Random number generators are designed to sample from specific distributions.

A random number generator that generates normally distributed random numbers is not identical to a random number generator that generates normally distributed random numbers (but excludes duplicate values)

Suppose you were using an RNG to generate random integers from a Poisson distribution with lambda = 10. Would you be happy if the RNG was “gimmicked” to exclude duplicate values? I’d argue that this would be a very flawed RNG because the samples being generated aren’t actually coming from a Poisson distribution.

Let’s move back to the example using “rand”. We’re dramatically changing the granularity of the problem, but the basic principle remains the same…

At the end of the day, I don’t think it’s reasonable evaluate an RNG based on a design criteria which is (arguably) contrary to the way in one would expect an RNG to behave (in this case, deliberately excluding duplicate values)

If you have an application that needs RNGs that are generated from “some distribution that looks very much like a normal distribution but is guaranteed never to produce a duplicate value” then you should probably find an RNG designed to do precisely that.

From my own perspective, I’d worry about implementations that were this persnicetky…

regards,

Richard

Walter Roberson
on 5 Feb 2011

Upper bound: 1/2+(1/2)*sqrt(1+8*S*ln(2)) where S is the number of distinct random numbers that can be generated.

Unfortunately Mathworks does not appear to document which Twister version they are using. Traditionally their random number generators were (0,1); mt19937ar over (0,1) is 32 bits precision, but over [0,1) is 53 bits precision.

If we speculate that Mathworks used the (0,1) version with 32 bits precision, then according to Wikipedia the 50% probability is at 7.7 × 10^4 numbers.

Peter Perkins
on 8 Feb 2011

Mark Shore
on 5 Feb 2011

Using the formula Walter gives leads to the following numbers for a 50% probability:

32-bit 7.716e4

53-bit 1.112e8

57-bit 4.470e8

64-bit 5.056e9

80-bit 1.294e12

A brute-force check of 1e10 random numbers

for n=1:1000 mindiff(n) = min(diff(sort(rand(1e7,1))))/eps(1); end

found 11 identical values, or a ~50% probability for 4.44e8 numbers, suggesting that the equivalent of a 57-bit precision RNG is in use.

Walter Roberson
on 6 Feb 2011

Walter Roberson
on 8 Feb 2011

The Matlab content of this question is:

- How many distinct random numbers can be generated in Matlab. (We still don't know for sure, but 2^53-1 appears likely.)
- Are the numbers generated with equal probability. (The algorithm would have to be examined in detail to be sure. If there is a bias, it is minute.)
- Is the period of the generator greater than the number of different possible values, or is it equal to the number of possible values?

Anything beyond that is a straight-forward application of the Birthday Paradox and is thus a mathematical question rather than a Matlab question. This forum exists to serve Matlab questions, not to serve mathematical questions and especially not to discuss the risks of fly-by-wire or to discuss deliberate fraud by Matlab users.

One only has to look at the period of MT and compare that to the largest number directly representable in Matlab in order to see that rand() will eventually generate a duplicate. Software that expects otherwise is broken.

It would be fair game to ask Mathworks to create a random number generator with no repetitions for those people who need one; it would even be fair game to ask Mathworks to rename rand() to make it clear that duplicates were possible. Those would be Matlab issues. Considering the amount of legacy code that uses rand(), I think Mathworks would want some strong evidence that "rand with replacement" was widely enough misunderstood to be worth renaming; I personally don't think speculation about possible trouble in genome software would qualify.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!