Probability 1/2 with MATLAB random number generators
조회 수: 17 (최근 30일)
이전 댓글 표시
Suppose you configure a MATLAB uniform random number stream (or use the default). Suppose you want two choices with equal probability, so 1/2 each. Without loss of generality, let the two outcomes be labeled 0 and 1.
Now, should the test be rand() > 0.5 or should the test be rand() >= 0.5 ?
mlfg6331_64, mrg32k3a, mt19937ar, dsfmt19937, and shr3cong are defined in terms of integral multiples of 2^(-N) for N = 53 or 52 or 32. There are 2^N possible N bit fractions that are exact multiples of 2^(-N), namely (0:2^N-1)/2^N . But 0 is excluded because it is the open interval, so the multiples are (1:2^N-1)/2^N . And there are always an odd number of those, just like with 1/10th's you have 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 -- 9 different possibilities, an odd number. If you include 0.5 exactly (a number which is precisely in the list of possibilities, being (2^(N-1))/(2^N) ) in the lower set then you give the lower set a bias; if you include 0.5 exactly in the upper set then you give the upper set a bias.
swb2712 indicates that it works natively in double precision and that all values in the open interval (0,1) are possible. If you typecast() non-negative double-precision numbers up to and including to 1.0 into uint64, you get the integers from 0 to typecast(1.0,'uint64') = 4607182418800017408 . That is an odd number of values including 0 and the last value corresponding to 1.0 exactly; when you exclude those two as being outside the open interval, you end up with an odd number of values and you are going to have bias when you divide it into two sets of equal probability. But notice that the half-way point through the integers is (4607182418800017408/2) = 2303591209400008704 so if you want to divide half way you should test for the value corresponding to that rather than testing for 0.5: the corresponding value is 1.11875110968003101149364479731944007560647073019006103414072774998094175776051774574042623356006417050422674048837541521801190921428379539907292472443881270372218911884781720662428061113996528944121698315475958721481213258908827606800475302345048341134463455488980253230058090875047749429690054193911503277704529859897485122299798376843682490289211273193359375e-154
mcg16807 is defined in terms of multiples of 1/(2^31 - 1) . The possibilities are (0:2^31-1)/(2^31-1) . The first of those gives 0, the last gives 1 exactly, so that would be (1:2^31-2)/(2^31-1) . There are an even number of these so you can divide into two sets with equal probability. But 0.5 exactly cannot occur in this scheme (not unless you try to push past 2^53, so whether you test for rand()>0.5 or rand()>=0.5 does not matter if you are using mcg16807 as your random number generator.
So... for all of the random number generators except swb2712 and mcg16807, testing for >0.5 or >=0.5 either way gives a bias and cannot reach probability 1/2 exactly. For swb2712 probability 1/2 exactly cannot be reached for the same reasons and the test should be near 1E-154 instead of 0.5 anyhow. mcg16807 alone can give probability 1/2 exactly, and for mcg16807 it does not matter if you use >0.5 or >=0.5 because it can not create exactly 0.5 .
Next question: how does randi() compensate for the number of generatable random values not being exactly divisible by the number of integers in the range to be generated?
댓글 수: 1
John Chilleri
2017년 1월 22일
The ever trivial and ever brainless unbiasing method:
x = rand();
while (x == .5)
x = rand();
end
As for the randi question. randi could be programmed similarly-brainlessly to eliminate bias:
% Assume n = 'number of generatable random values'
% Assume k = 'number of integers in the range'
% Let rv be the randomly generated value
UBound = floor(n/k)*k;
while (rv > Ubound)
Generate new rv;
end
% Proceed
If you're confused as to why I'm randomly commenting on this, I just wanted to browse the questions of number 1.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Random Number Generation에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!