bugGNU Octave - Bugs: bug #54619, randi() is biased

 
 

bug #54619: randi() is biased

Submitter:  None
Submitted:  Tue 04 Sep 2018 04:40:58 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Cris Luengo Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 27 Mar 2023 09:37:20 PM UTC, comment #29: 

Thank you for checking.  Replacement of 'flintmax + 1' with 'flintmax' done in this changeset http://hg.savannah.gnu.org/hgweb/octave/rev/9d5ccdbae3ad.

Rik <rik5>
Group administrator
Mon 27 Mar 2023 09:26:53 PM UTC, comment #28: 

I cannot remember why I thought this flintmax()-1 was needed, and now I am actually quite convinced that it gives incorrect results, as the line


r_prim = floor ( rand (M,1) * (flintmax ()-1) );


would lead to 0 appearing in r_prim with double probability. No, using just flintmax() everywhere instead of flintmax()+1 should be perfectly okay -- the results should stay perfectly the same, but the reader of the code would not be confused.

By the way, I will document what I found about octave's rand() to answer my questions in comment #19: according to https://hal.science/hal-02427338, octave's rand() works by first generating a uniformly distributed integer between 0 and 2^53-1, adds 0.4 to it (looks like an ad hoc-decision, as no other popular software package does it that way) and divides by 2^53. Thus, it obviates (partly?) the problem that the less-significant bits of the returned doubles below 2^-11 are all zero, at the price of a very small bias towards lower values.

In the context of randi(), everything should be perfect: multiplication by flintmax() undoes the division by 2^53, floor undoes the addition of 0.4, so that r_prim should really be a uniformly distributed integer. Thus, apart from the potential issue of inexact floating-point division as mentioned in comment #19, as far as I can tell randi() should work as advertised.

Michael Leitner <mleitner>
Mon 20 Mar 2023 12:09:08 AM UTC, comment #27: 

@mleitner: Someone has asked over on the Octave Discourse site (https://octave.discourse.group/t/baffling-expressions-in-randi-implementation/4138) why the code in randi.m has expressions of this form:


(flintmax() + 1)


It is true that flintmax() is the largest representable integer and that


(flintmax() + 1) == flintmax()


so it would seem that there is a simplification that could be made.  But, I wanted to ask here since you wrote the code and might have had a reason that I don't immediately divine.

Rik <rik5>
Group administrator
Wed 19 Dec 2018 06:06:19 AM UTC, comment #26: 

I reviewed and amended Michael's patch and pushed it here (https://hg.savannah.gnu.org/hgweb/octave/rev/25d3e8e49d5c).  This will be a part of the 5.0 release in early January.  Marking as fixed and closing report.

Rik <rik5>
Group administrator
Mon 10 Sep 2018 11:17:17 PM UTC, comment #25: 

@Michael: All fair comments.  I saw that the loop would be executed infrequently, but didn't continue to calculate how infrequently that might be.

The worst case is unbounded.  But, what I had in mind was column 3 of Table 1, "Maximal number of remainders per integer" which is Infinity for the Java algorithm which is why I wasn't interested in it.

Total time = time_to_create_random_integer + time_map_integer_onto_range. 

Running the openbsd algorithm in randi_rej.cpp with


for i = 1:20
tic;
x=randi_rej (6004799503160661);
bm(i)=toc;
end

octave:19> mean (bm(2:end))
ans =  0.026421


If I change the code to


  for (octave_idx_type i = 0; i < sz; i++)
    {
      data[i] = octave::randi53 ();
      //data[i] = openbsd (s) + imin;
    }


, re-compile, and re-run the 20 point benchmark I get


octave:3> mean (bm(2:end))
ans =  0.011396


So the random number generation is about 43% of the total time and the mapping is 57% of the time.  Since it is so balanced, I think there are opportunities for improvement in both halves.  If I run the nearlydivisionless algorithm on a range where it works, imax=9, the results are impressive.


octave:7> mean (bm(2:end))
ans =  0.018590


But, since it didn't work for all values of imax I can't advocate for it.

Rik <rik5>
Group administrator
Mon 10 Sep 2018 09:51:20 PM UTC, comment #24: 

Rik, I have to comment on a number of points in your series of comments.

First, I would say that there are just two main variants of a single algorithm for deriving unbiased random integers in an arbitrary range from a source of unbiased random integers in a given range. Both use rejection: For instance, when you have a generator for integers in [0,65535] (i.e., 16 bit) and you want to generate integers in [0,9], you first generate a source integer, which you reject and regenerate iteratively if it should fall within [65530,65535]. So now you have an unbiased integer r within [0,65529], and you get your output by doing either mod(r,10) or floor(r/6553) (when you are doing integer arithmetic, you can drop the floor). What is discussed in the arxiv-manuscript is how you can do that while keeping the number of integer divisions low. So actually all the implementations have unbounded worst-case behaviour, because they all use rejection.

I think you misunderstood my solution. Yes, there is a loop, but it does not loop over the elements. Instead, I generate first more numbers that I expect to need, do everything vectorized, and only when at the end I see that the number of generated numbers is not large enough because I had to reject too many, I try again. If have the comment "should practically always be true" in the code, because I generate use an excess number that is 10 times the expected standard deviation of the number of rejections. In octave, normcdf(10) evaluates to exactly 1, while 1-normcdf(8) is 1e-15. Thus the probability for the loop having to be evaluated a second time is probably about 1e-20. That's what I mean by "practically always".

Yes, a compiled function will always be faster than a vectorized m-file, mainly because of the fact that with vectorization, you have to repeatedly retrieve and store back your intermediate vectors to main memory, while in a loop you have only to store the final value, everything else happens in registers. Further, when you directly use randi53(), you obviate the cast to double and subsequent division that happen in rand(), and you do not have to undo the division. I think that these are the main reasons why your code is faster, not so much that the arxiv-algorithms are better (because they are actually the same).

So how fast is my code? I see


> tic;a = randi(6004799503160661,1000000,1);toc
Elapsed time is 0.245796 seconds.
> tic;a = randi_(6004799503160661,1000000,1);toc
Elapsed time is 0.547995 seconds.
> tic;a = randi_(10,1000000,1);toc
Elapsed time is 0.398885 seconds.


About a factor of 2.2 slower than the simple biased m-file when one in three values has to be rejected, dropping to a factor 1.6 when practically nothing is rejected. The slowing-down with rejection cannot be avoided, and a factor of 1.6 for an m-file is not bad, I think. So as long as it is about addressing the bias, and as long as nobody finds a bug in my implementation, I still think it can be used.

Thinking about what to do when you really want to beat Matlab: I think you meant that rand should keep returning doubles within [0,1). Yes, that goes without saying. For this you need randi53. But what would be nice to expose would be something that indeed returns a uint64 with the full range of [0,2^64-1], because this is in most likelihood what the Mersenne twister internally returns. And further, for beating Matlab with something "new" I would propose to change the underlying PRNG to something of the PCG family http://www.pcg-random.org/, which are both among the fastest and "randomest" PRNGs available today. I see that Melissa O'Neill even has a very recent post about just the problem we are discussing here: http://www.pcg-random.org/posts/bounded-rands.html


Michael Leitner <mleitner>
Mon 10 Sep 2018 06:00:30 PM UTC, comment #23: 

Going forward, I have other bug reports that I want to sqaush so I won't be working on this.  My suggestion would be implement the functionality in liboctave in oct-rand.cc or a combination of oct-rand.cc and randmtzig.cc.  The interpreter then needs a way to access the new functionality.  I would probably put that in libinterp/corefcn/rand.cc.  There is code in rand.cc, or fill_matrix() in data.cc, for extracting the dimensions from a variable number of input arguments which could be borrowed for randi input processing.

For the moment I would prefer to keep rand returning values within [0,2^53).  The base data type in Octave is a double with only 53 bits of precision and I fear introducing more bias during the conversion from 64 to 53 bits.  Also, minimal change necessary to solve the problem seems best from a software engineering perspective.

There is still a good GSOC project here to convert to using C++ standard libraries for random numbers.  Maybe simultaneously with that project we could investigate the impact of returning a uint64_t from an internal API.

Rik <rik5>
Group administrator
Mon 10 Sep 2018 05:31:00 PM UTC, comment #22: 

It would be fun to have bragging rights to an algorithm which is both faster than Matlab and unbiased so I tried implementing the "new" algorithm in the paper in randi_new.cpp.  The author outlined the algorithm, and also provided a code implementation for L=64.  Because we have L=53 I could not use his code directly and instead coded the algorithm as I understood it.  What I found was that it worked for small values of imax, but hung in an essentially infinite loop for large values.  That may be my misunderstanding, but I wasn't super motivated to find out.  The algorithm relies on 128-bit unsigned data types, which g++ does support, but they are not part of the C++ standard and so there would need to be some extra work in the configure script to determine which compiler was being used, whether 128-bit primitives were supported, and a fallback mechanism if they were not.

Rik <rik5>
Group administrator
Mon 10 Sep 2018 05:01:10 PM UTC, comment #21: 

For testing, I used a script tst_randi.m which is shown below, and attached.


#a = randi (6004799503160661, 1e6, 1);
#a = randi (9,1e6,1);
#a = randi_rej (6004799503160661);
#a = randi_rej (9);
#a = randi_new (6004799503160661);
a = randi_new (9);

b = a(a>median(a));
c = hist(mod(b,2),[0 1])
c ./ sum (c)


It generates one million values using a particular method and then uses Michael's test even/odd test.  It prints out both the absolute counts and the percentages (which should be 50%:50%).

As a warmup, running tst_randi with the existing algorithm (randi) produces


octave:13> tst_randi
c =

   166430   333570

ans =

   0.33286   0.66714


So, there is a clear bias in the results using the existing algorithm.  For some performance benchmarkings,


octave:14> tic; x = randi (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0255809 seconds.
octave:15> tic; x = randi (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0356832 seconds.
octave:16> tic; x = randi (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0255952 seconds.
octave:17> tic; x = randi (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0247469 seconds.
octave:18> tic; x = randi (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0325689 seconds.


So, roughly, 30 milliseconds.  Using the rejection method,


octave:2> tic; x = randi_rej (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0560811 seconds.
octave:3> tic; x = randi_rej (6004799503160661, 1e6, 1); toc
Elapsed time is 0.0398862 seconds.
octave:4> tic; x = randi_rej (6004799503160661, 1e6, 1); toc
Elapsed time is 0.040313 seconds.
octave:5> tic; x = randi_rej (6004799503160661, 1e6, 1); toc
Elapsed time is 0.040343 seconds.
octave:6> tic; x = randi_rej (6004799503160661, 1e6, 1); toc
Elapsed time is 0.046854 seconds.


Very roughly, 45 milliseconds, or about 50% slower.  Quoting from the paper,


The number of random integers consumed by this algorithm follows a geometric distribution, with a success probability p = 1 − (2^L mod s)/2^L . On average, we need 1/p random words. This average is less than two, irrespective of the value of s. The algorithm always requires the computation of two remainders.


The worst case is when 2^L mod s is large which is at (2^L / 2) + 1.  For our case, this is (2^53 / 2) + 1 = 4503599627370497.  Testing this idea proves to be correct.


octave:2> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.0646839 seconds.
octave:3> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.102659 seconds.
octave:4> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.0820141 seconds.
octave:5> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.0970309 seconds.
octave:6> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.0697341 seconds.
octave:7> tic; x = randi_rej (4503599627370497); toc
Elapsed time is 0.0705478 seconds.


Maybe this averages 80 milliseconds.  But of course, the real reason to use thismethod is because it is unbiased.  Running tst_randi shows


octave:8> tst_randi
c =

   250205   249795

ans =

   0.50041   0.49959


Also, it is unlikely that users are choosing pathological values for imax frequently.



(file #44968)

Rik <rik5>
Group administrator
Mon 10 Sep 2018 04:33:15 PM UTC, comment #20: 

I did some work on this problem over the weekend just to hack something up.  First, I'm attaching the academic paper that Michael Godfrey pointed to (Fast_Random_Integer_Generation_in_an_Interval.pdf).  This outlines two methods used by existing libraries: 1) rejection (OpenBSD, GNU libc, others), 2) Java.  Of these, rejection is used by many and it has guaranteed worst case behavior.  The Java library has better average performance, but unbounded worst case behavior, which I'm not interested in adopting.  Finally, the author proposes a third method of his own devising which I will call "new" method.  The paper is from May, 2018 so this is very "new" compared to the existing methods.

I appreciate Michael Leitner's m-file approach, but it does contain a do-until loop.  Because Octave is interpreted, and we don't have a JIT implementation for loops, any looping behavior involves a big performance hit.  In situations where looping can't be avoided the suggestion is to move to C++ which is what I have done.

Second attachment is a small change to expose the function randi53 (53-bit random integer) from liboctave/numeric/randmtzig.cc.  See randi53.diff.  Because of some vagaries of the build system to use this you will need to patch the Octave source tree and then install Octave so that 'mkoctfile' has the new header file randmtzig.h.  Instructions


cd octave_src_dir
patch -p1 < randi53.diff
make
make install


After that you can work with the two octfile samples I made: randi_rej.cpp and randi_new.cpp.  These are not complete--there is no input validation, documentation, support for specifying number of values, BIST tests, etc.--but they do allow you to set imin and imax and they return a column vector [1e6,1] of values.

To compile them use


mkoctfile -v -O2 randi_rej.cpp
mkoctfile -v -O2 randi_new.cpp


More to follow, but I can only attach 4 files at a time so I have to continue in a new comment.



(file #44964, file #44965, file #44966, file #44967)

Rik <rik5>
Group administrator
Sun 09 Sep 2018 09:49:29 PM UTC, comment #19: 

See the attached patch (a diff against the current randi.m). I think it is correct, at least it passes the test on the ratio of evens and odds. This could be also added in the tests, but I do not know how to write a test on stochastic functions. Of course, that it does pass this specific test does not prove anything, so I invite also other people to have a look on it. Essentially it just tries to create integers between 0 and flintmax-1 (2^53-1 for IEEE-754 doubles), divides them by another integer K so that the requested range is covered (for exact integer arithmetic this should result in exactly uniform probabilities within the requested range) and rejects all above the range. It initially generates a buffer of random numbers that is in all likelihood larger than what is requested, so that it is only insignificantly slower than the previous implementation for cases where not many numbers have to be dropped (for small ranges, that is, where the previous implementation was already nearly correct). In these cases, the sequence of returned numbers is also the same, provided not too many numbers are requested.

A potential issue: is IEEE-754 double division predicted to be exact on exactly representable integers, when the result is again an integer? If division is multiplication by the inverse, then probably not, because the inverse is in general not exactly representable. In this case, it would probably be better to do these manipulations in uint64, which should be required to be exact.

As a side note: I would be very much in favour of exposing the fundamental random number generator not via rand(), which returns doubles in the range [0,1], but by a function returning uniformly distributed uint64, which they internally indeed are. This would save operations in the present problem which have only the purpose to undo this conversion to scaled doubles, further you can get the scaled doubles by a two-line m-file, and finally this would then document what rand() really returns (is exact 0 and/or exact 1 included, if so, with which probability, and is indeed every double between 0 and 1 returned or only those that are rounded uint64 divided by 2^64 and so on).

(file #44961)

Michael Leitner <mleitner>
Sat 08 Sep 2018 06:39:52 PM UTC, comment #18: 

Thanks Michael for the observable example.  I've changed the release field to dev.  This should be fixed for the 5.0 release due at the end of this year.

Rik <rik5>
Group administrator
Sat 08 Sep 2018 06:28:24 PM UTC, comment #17: 

dev Octave (~4 sept):

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


...and Matlab r2018b prerelease:

>> a = randi(6004799503160661,10000,1,'uint64');
Error using randi
Outputs of class uint64 not supported.


Philip Nienhuis <philipnienhuis>
Group Member
Sat 08 Sep 2018 04:09:36 PM UTC, comment #16: 

I think the test case should be changed to:


a = randi(6004799503160661,10000,1,'uint64');


Thanks, rahnema1

Anonymous
Sat 08 Sep 2018 08:55:02 AM UTC, comment #15: 

So at least it's a Matlab-compatible bug

Philip Nienhuis <philipnienhuis>
Group Member
Sat 08 Sep 2018 08:53:42 AM UTC, comment #14: 

Matlab:

>> version
ans =
    '9.5.0.882065 (R2018b) Prerelease'

>> a = randi (6004799503160661, 10000, 1);

>> b = a(a>median (a));

>> c = hist (mod (b,2), [0 1])
c =
        1640        3360


Philip Nienhuis <philipnienhuis>
Group Member
Sat 08 Sep 2018 01:36:09 AM UTC, comment #13: 

Thanks for illuminating!

Octave

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

   1689   3311


Julia seems to get it right

using StatsBase
a = rand(1:6004799503160661, Int(1e4));
b = a[a .> median(a)];
h = fit(Histogram, mod.(b,2), nbins=2, closed=:right)

Histogram{Int64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  -0.5:0.5:1.0
weights: [2474, 0, 2526]
closed: right
isdensity: false

Note: there seems to be abug in fit which fits 3 bins instead of 2

Numpy aslo gets it right

import numpy as np
a = np.random.randint(low=1,high=6004799503160661, size=int(1e4))
b = a[a > np.median(a)]
nn, xx = np.histogram(np.mod(b,2), 2)
print(nn, xx)
[2524 2476] [0.  0.5 1. ]



From my view this is now a serious bug, we should fix it asap.

Juan Pablo Carbajal <juanpi>
Group Member
Sat 08 Sep 2018 12:01:25 AM UTC, comment #12: 

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.

Michael Leitner <mleitner>
Fri 07 Sep 2018 04:33:16 PM UTC, comment #11: 

The code below to try to quantify the bias is pointless, sorry. Verifying correctness of a PRNG is a very complicated thing and nobody knows how to do it right. You need to look at this issue from first principles: compute the probability for each output value.

Here's a quick explanation of the problem with a extreme example:

Assume you have a 2-bit RNG (generates values in the range [0,3]. You want to generate integers in the range [1,3]. You use Octave's method:


double x = RNG / 4
    # x = 0.00 -- 25%
    # x = 0.25 -- 25%
    # x = 0.50 -- 25%
    # x = 0.75 -- 25%
int n = floor(x * 3) + 1
    # x = 0.00 -> n = 1
    # x = 0.25 -> n = 1
    # x = 0.50 -> n = 2
    # x = 0.75 -> n = 3


You are now mapping 4 possible values into 3 bins. 2 of these bins will have a 25% chance of being hit, 1 bin will have a 50% chance. You can do the mapping however else you want, the outcome will always have this bias. The only solution is to reject 1 in 4 results of the RNG.

With a 53-bit RNG as Octave uses, this effect is much, much smaller, but it is still there.

I hope this example clarifies the problem.

I suggested earlier to "throwing away randmtzig.cc and replacing that with stuff from the Standard Library", but what I really meant was to replace the `randi` M-file with a call to the Standard Library. The random number generator you use is OK, and it might even be possible to use `std::uniform_int_distribution` with it. I would recommend replacing that also, just because it's less code to maintain, but if you need to keep the old behaviour intact, maybe make the fewest changes that fix `randi`.

Anonymous
Fri 07 Sep 2018 04:14:03 PM UTC, comment #10: 

Rik,

Sounds good to me. This is a good Summer project task for a start.
The C++ code should in any case be selected based on an option
setting. As you said, the old code will need to be supported for
compatibility reasons. That also means that the sooner the
old code is frozen the better. This also supports (my) view
that "fixing" the bias in the current code is not a good idea.

Michael Godfrey <godfrey>
Group Member
Fri 07 Sep 2018 03:27:09 PM UTC, comment #9: 

@Michael: I think this might make a good Google Summer of Code project.  This task isn't large enough to be a standalone project, but it could be part of a larger project to change over the core to use the C++ libraries for random number generation.  I read through the paper you sent and I was ambivalent.  It's quite nicely worked out, and there are code samples, but if we move in that direction then we also are back in the business of maintaining a code base for random numbers.  There is a simplicity, and reduced maintenance burden, to just calling library code and letting the library maintainer's worry about such things.  Of course, a true engineering decision would look at a number of variables like accuracy, performance, maintenance burden, etc.  That's why I think it might make a good GSOC project for someone to explore these concerns and then implement a solution after consultation with the Octave Maintainers.

Rik <rik5>
Group administrator
Fri 07 Sep 2018 09:27:05 AM UTC, comment #8: 

Rik,

Discussion of "correctness" of random number generators
has a long history and has not been particularly helpful.
The idea of a "perfect" random number generator appears
to be wishful thinking.
However, using the C++ implementation for this case seems
a reasonable choice combined with the option of backward
compatibility as you point out.

Michael Godfrey <godfrey>
Group Member
Thu 06 Sep 2018 09:39:31 PM UTC, comment #7: 

Rik,

This deserves a careful look. The paper at:
https://arxiv.org/pdf/1805.10941.pdf

may also be helpful.

Michael Godfrey <godfrey>
Group Member
Thu 06 Sep 2018 09:00:17 PM UTC, comment #6: 

For the sake of comparisson I run Rik's code in Julia 1.0.0, I get thefollowing deviations


 0.01408
 0.04265
 0.0012799999999999999
 0.01442
 0.00066
 0.00251
 0.01893
 0.00568
 0.0005200000000000001
 0.020790000000000003

max 0.04% in bin 2

and a second time

 0.00303
 0.03839
 0.04325
 0.02103
 0.05595
 0.028669999999999998
 0.03591
 0.025109999999999997
 0.00456
 0.0053

max 0.056% in bin 5

In numpy 1.14.2 I get

array([0.00682, 0.0383 , 0.02544, 0.02562, 0.05767, 0.01353, 0.07619, 0.04365, 0.00261, 0.01745])

max 0.076% at bin 7
and again

array([0.00443, 0.05387, 0.01102, 0.06105, 0.0012 , 0.00123, 0.04119, 0.02084, 0.07958, 0.02007])

max 0.08% at bin 9

Both are similar to Octave's deviations.

I appreciate the comments from the OP. If he can just illustrate the value of his suggestion with a working example, then I am pretty somebody will be happy to improve the method.
It seems there is no agreement on what "serious scientific software" is, at least Julia and numpy are considered pretty serious.

Julia code:

using StatsBase
x = rand(1:10, Int(1e8));
h = fit(Histogram, x, nbins=10, closed=:left)
errpct = abs.(h.weights .- 1e7) / 1e7 * 100


Numpy code:

import numpy as np
x = np.random.randint(low=1,high=11, size=int(1e8))
nn, xx = np.histogram(x)
errpct = np.abs(nn - 1e7) / 1e7 * 100


Juan Pablo Carbajal <juanpi>
Group Member
Thu 06 Sep 2018 08:57:46 PM UTC, comment #5: 

I hear the assertion that the Octave method is incorrect, but I haven't been able to determine that myself.  But, assuming it is incorrect, is it inadequate?  That's why I want a numerical estimate of the bias of the current system.  The last argument, not re-inventing the wheel, is the most compelling.  Octave only recently began requiring C++11 so it is now possible to contemplate using the STL for this case.  Long term, we would like to move all of the random number generation to the standard library, but we also have to support legacy random number generators (as Matlab does), so the code is always going to be rather complicated so that programmers can reproduce simulations they have written years ago.

Rik <rik5>
Group administrator
Thu 06 Sep 2018 08:26:03 PM UTC, comment #4: 

There is a correct way to generate random integers in a given range, and there are incorrect ways. The method that I linked is the correct way. Other ways are incorrect. The bias might be small, but it's there.

GNU's C++ Standard Library does it right also:


const __uctype __uerange = __urange + 1; // __urange can be zero
const __uctype __scaling = __urngrange / __uerange;
const __uctype __past = __uerange * __scaling;
do
   __ret = __uctype(__urng()) - __urngmin;
while (__ret >= __past);
__ret /= __scaling;


No serious scientific software would use the method used by Octave. This is why I'm reporting it as a bug. Take it or leave it.

Note that random floating-point numbers are totally different from random integers, and should not be used as a comparison. And they should not be used as a basis either.

This "large project" is just throwing away randmtzig.cc and replacing that with stuff from the Standard Library. Random number generation is complicated enough that one should not attempt to reinvent the wheel.

Anonymous
Thu 06 Sep 2018 06:24:07 PM UTC, comment #3: 

the rand from Octave is implemented in liboctave/numeric/randmtzig.cc.  As far as I can tell, there is no 64-bit integer that is cast to double.  Similarly, there is no "division + floor" in randi.  There is multiplication, and floor is used to remove the fractional part.  Because I can't verify some of the first assertions, I can't go on to accept the conclusions that follow.  The Stack Overflow response has this comment


The suggestion to use floating-point division is mathematically plausible but suffers from rounding issues in principle. Perhaps double is high-enough precision to make it work; perhaps not. I don't know and I don't want to have to figure it out; in any case, the answer is system-dependent.


So it's not clear, even to the original author, that Octave's current system doesn't work well enough.  Before any large project is undertaken it makes sense to have an idea of the costs and the benefits.  It seems to me that the benefits need to be quantified.  How large is the bias of the existing system (if it even exists) that is implemented in Octave?

Also, relevant, are we any worse than the competition?  According to The Mathworks, the algorithm Octave implements is the one suggested for random double numbers (http://www.mathworks.com/help/matlab/ref/rand.html)


Generate a 10-by-1 column vector of uniformly distributed numbers in the interval (-5,5).

r = -5 + (5+5)*rand(10,1)

r = 10×1

    3.1472
    4.0579
   -3.7301
    4.1338
    1.3236
   -4.0246
   -2.2150
    0.4688
    4.5751
    4.6489

In general, you can generate N random numbers in the interval (a,b) with the formula r = a + (b-a).*rand(N,1).


For random integers, The Mathworks suggestion is to use randi(), but I don't know if they implement something special there or not.

As an attempt at quantifying the bias, I used


x = randi ([1, 10], 1e8, 1);
[nn,xx] = hist (x, 1:10);
errpct = abs (nn - 1e7) / 1e7 * 100
errpct =

 Columns 1 through 8:

   0.0376100   0.0484100   0.0173600   0.0174000   0.0120900   0.0249200   0.0092400   0.0054200

 Columns 9 and 10:

   0.0250300   0.0050400


So, as much as .05% deviation, but this could be statistical fluctuations.  Running a second time I get


errpct =

 Columns 1 through 8:

   0.0057600   0.0249700   0.0160600   0.0201400   0.0295600   0.0249300   0.0148500   0.0361300

 Columns 9 and 10:

   0.0643700   0.0242500


This time it is a .06% deviation, but it is not in the same bin as before so it doesn't look like systematic bias.


Rik <rik5>
Group administrator
Thu 06 Sep 2018 04:15:06 PM UTC, comment #2: 

This has nothing to do with the C `rand`. Octave's `randi` is an M-file, you can type `edit randi` to see it. It uses Octave's `rand`. Line 106:


ri = imin + floor ( (imax-imin+1)*rand (varargin{:}) );


`rand` uses the Mersenne Twister algorithm, which outputs a 64-bit   integer, cast to a double. Even if this cast is done correctly, leading to a perfectly uniform distribution of values, this double float value doesn't sample integer values equidistantly, and the output range likely doesn't divide the input range evenly either. Thus, the division + floor in `randi` will lead to a bias in the output distribution.

In short, a correct `randi` is implemented using the integer output of the Mersenne Twister, and using the rejection scheme as shown in the linked Stack Overflow answer.

Anonymous
Wed 05 Sep 2018 10:52:41 PM UTC, comment #1: 

Could you point to the code in Octave where this division occurs?  Octave is not using rand() or random() from the stdlib for its implementation of rand() at the interpreter level.  Thus, I'm not sure the analysis applies.

Rik <rik5>
Group administrator
Tue 04 Sep 2018 04:40:58 PM UTC, original submission:  

The randi function uses an incorrect algorithm that leads to biased results. It is implemented in terms of `rand`. But `rand` is created by dividing the output of a 64-bit random number generator by 2^64-1, scaling this value up and rounding down doesn't lead to a uniform distribution, as the original integer intervals are most likely not evenly divisible into the output interval, and thus some output values are slightly more likely than other output values.

This is maybe only a small bias, but it is there.

See here for a correct algorithm: https://stackoverflow.com/questions/2509679/how-to-generate-a-random-integer-number-from-within-a-range/6852396#6852396

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #44968:  tst_randi.m added by rik5 (229B - text/x-matlab)
file #44965:  randi_rej.cpp added by rik5 (841B - text/x-c++src)
file #44966:  randi53.diff added by rik5 (1KiB - text/x-patch)
file #44967:  randi_new.cpp added by rik5 (1KiB - text/x-c++src)
file #44961:  bug54619_patch added by mleitner (966B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-12-19 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2018-09-10 rik5 Attached File- Added tst_randi.m, #44968
    2018-09-10 rik5 Attached File- Added Fast_Random_Integer_Generation_in_an_Interval.pdf, #44964
        Attached File- Added randi_rej.cpp, #44965
        Attached File- Added randi53.diff, #44966
        Attached File- Added randi_new.cpp, #44967
    2018-09-09 mleitner Attached File- Added bug54619_patch, #44961
    2018-09-08 rik5 StatusNeed Info Confirmed
        Release4.0.0 dev
    2018-09-05 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code