From r@lf@@tubner @ending from d@q@n@@com  Thu Sep 20 12:59:00 2018
From: r@lf@@tubner @ending from d@q@n@@com (Ralf Stubner)
Date: Thu, 20 Sep 2018 12:59:00 +0200
Subject: [Rd] Bias in R's random integers?
In-Reply-To: <CAN_1p9znLss3QKuF3WB5ucM_F2=paNx28eCdJecB8fj0p_w-6Q@mail.gmail.com>
References: <CAN_1p9wVX=AR8aYtC2LJo0bgN9cRW3XcbYJqnVHdCzA-QA1ceg@mail.gmail.com>
 <ffb74021-8d95-a820-8013-f0e51dc1a9b3@gmail.com>
 <CALEXWq3booRRTSYvYqcjV7qPZGOC1oh5o45b5o=J2ypQi5s-ag@mail.gmail.com>
 <ece52410-ddb2-2198-9291-201475d3f46e@gmail.com>
 <CAA7KnHHEejYVn8NtKrY1KuKcJUq8n4LoJSDx9BZv-LJyC=r9yQ@mail.gmail.com>
 <73c4667d-38ad-f5ef-468c-4c368e0fdc12@gmail.com>
 <CAA7KnHFcY5YJ2kTAfPJJ0PHwcWx7CumE49A4NgJWC2NK5CV=-A@mail.gmail.com>
 <c3590da8-2e92-4e46-3a81-0d0dbfc81845@gmail.com>
 <CAA7KnHE9OSxML3Cp8Y-gd7CbBj-86UP9qAy6yzgvfcb+=Vg86Q@mail.gmail.com>
 <e82cfff0-3a8e-1d7b-de5c-02b69974ae16@gmail.com>
 <4143568b-c67f-2309-f9ec-676f9e791747@gmail.com>
 <CAN_1p9znLss3QKuF3WB5ucM_F2=paNx28eCdJecB8fj0p_w-6Q@mail.gmail.com>
Message-ID: <5245f770-a010-b976-28c8-13a2ec74a088@daqana.com>

On 9/20/18 1:43 AM, Carl Boettiger wrote:
> For a well-tested C algorithm, based on my reading of Lemire, the unbiased
> "algorithm 3" in https://arxiv.org/abs/1805.10941 is part already of the C
> standard library in OpenBSD and macOS (as arc4random_uniform), and in the
> GNU standard library.  Lemire also provides C++ code in the appendix of his
> piece for both this and the faster "nearly divisionless" algorithm.
> 
> It would be excellent if any R core members were interested in considering
> bindings to these algorithms as a patch, or might express expectations for
> how that patch would have to operate (e.g. re Duncan's comment about
> non-integer arguments to sample size).  Otherwise, an R package binding
> seems like a good starting point, but I'm not the right volunteer.
It is difficult to do this in a package, since R does not provide access
to the random bits generated by the RNG. Only a float in (0,1) is
available via unif_rand(). However, if one is willing to use an external
RNG, it is of course possible. After reading about Lemire's work [1], I
had planned to integrate such an unbiased sampling scheme into the dqrng
package, which I have now started. [2]

Using Duncan's example, the results look much better:

> library(dqrng)
> m <- (2/5)*2^32
> y <- dqsample(m, 1000000, replace = TRUE)
> table(y %% 2)

     0      1
500252 499748

Currently I am taking the other interpretation of "truncated":

> table(dqsample(2.5, 1000000, replace = TRUE))

     0      1
499894 500106

I will adjust this to whatever is decided for base R.


However, there is currently neither long vector nor weighted sampling
support. And the performance without replacement is quite bad compared
to R's algorithm with hashing.

cheerio
ralf

[1] via http://www.pcg-random.org/posts/bounded-rands.html
[2] https://github.com/daqana/dqrng/tree/feature/sample

-- 
Ralf Stubner
Senior Software Engineer / Trainer

daqana GmbH
Dortustra?e 48
14467 Potsdam

T: +49 331 23 61 93 11
F: +49 331 23 61 93 90
M: +49 162 20 91 196
Mail: ralf.stubner at daqana.com

Sitz: Potsdam
Register: AG Potsdam HRB 27966 P
Ust.-IdNr.: DE300072622
Gesch?ftsf?hrer: Prof. Dr. Dr. Karl-Kuno Kunze


-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 833 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20180920/cce85e63/attachment.sig>


[Rd] Bias in R's random integers?
Ralf Stubner