Does MATLAB have a Birthday Problem?
조회 수: 9 (최근 30일)
이전 댓글 표시
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.
댓글 수: 5
Malcolm Lidierth
2011년 2월 8일
In the UK, when the BBC first moved its football highlights show ('Match of the Day') from its traditional Saturday night slot, the birth rate tripled 9 months later. The birthday paradox, makes a false assumption: birthdays are not random - they vary seasonally and are influenced by all sorts of things.
채택된 답변
Derek O'Connor
2011년 2월 5일
댓글 수: 10
Walter Roberson
2011년 2월 6일
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.
추가 답변 (4개)
Richard Willey
2011년 2월 7일
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
댓글 수: 4
Richard Willey
2011년 2월 9일
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
2011년 2월 5일
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.
댓글 수: 1
Peter Perkins
2011년 2월 8일
Walter, the User Guide includes this: "mt19937ar: The Mersenne Twister, as described in [9], with Mersenne prime 2^19937-1. This is the generator documented at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html. It has a period of 2^19937-1. Each U(0,1) value is created using two 32-bit integers from the generator; the possible values are all multiples of 2^-53 strictly within the interval (0,1)."
Mark Shore
2011년 2월 5일
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.
댓글 수: 14
Bruno Luong
2011년 2월 6일
@Jan: Walter is right by stating the fact that Matlab rand() command returns only integer-multiple of 2^(-53).
Walter Roberson
2011년 2월 6일
Jan, the testing I did was unable to find a non-zero difference smaller than 2^(-53) as well. But it could just be quite rare. I don't think I have time to generate 2^53-ish values to check!
Walter Roberson
2011년 2월 8일
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.
댓글 수: 10
Peter Perkins
2011년 2월 10일
While the code on the FEX is indeed a wrapper around N&M's original mt19937ar code, the wrapper differs slightly from the implementation that the RAND function has used for the MT since that generator was officially supported circa R2007b. You have noticed two things: the range, and the speed. There is also the treatment of a zero seed. I can't recall if there are other differences, but in any case, RAND is RAND, and the FEX submission is the FEX submission.
The doc for RAND and my earlier comments are correct, and in particular RAND (internally) rejects an exact (double precision) zero in the very infrequent event that it happens to generate one.
Hoping this clears it all up. I almost hesitate to even mention it, but as described in the RandStream doc, you can also set the FullPrecision property of an mt19937ar stream to false, and get values that are multiples of 2^-32.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!