[Octave-bug-tracker] [bug #54619] randi() is biased

[Michael Leitner]

Follow-up Comment #12, bug #54619 (project octave):

Under normal operation, the bias is really tiny, so that you will never be
able to see it with such tests. The deviations you see are just what you
expect for a binomial distribution. But as we know how the code is working, we
can easily construct a pathological problem:


a=randi(6004799503160661,10000,1);
b=a(a>median(a));
c=hist(mod(b,2),[0 1])


This just generates a bunch of random integers, takes the larger half of them,
and computes how many of them are even or odd. The result should be that on
average exactly half of them are odd, but here you see that it is rather two
thirds. So this is a significant deviation, and when you blindly assume that
randi is doing its job correctly, you could generate capital nonsense. 

This works because the magic number 6004799503160661 is bitmax()*2/3, so when
we assume that rand() itself works correctly, that is, that every possible
double between 0 and 1 appears with the correct probability (not with equal
probability, because there is exactly the same number of floating-point
numbers between 0.25 and 0.5 as between 0.5 and 1), then in the upper half of
numbers either one or two doubles get assigned to one and the same integer.

It would be interesting to see how Matlab performs here, or Julia and numpy.
If most of the competition does it correctly, then I would just quickly hack
up the rejection part, which would not be too difficult.

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?54619>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/