Vor etwa einem halben Jahr beschrieb ich dqsample als eine faire Alternative zu base::sample(). Inzwischen ist das nur mehr von historischem Interesse, nachdem die R Version 3.6.0 einen aktualisierten Algorithmus für Zufallsstichproben enthalten wird. Einige der Techniken aus dqsample sind nun aber Teil von dqrng. Selbstverständlich ist der verwendete Algorithmus weiterhin fair:
|
1 2 3 4 5 |
m <- 2/5 * 2^32 x <- dqrng::dqsample(floor(m), 1000000, replace = TRUE) plot(density(x[x %% 2 == 0]), main = "dqrng::dqsample", xlab = NA) lines(density(x[x %% 2 == 1]), col = "#FF8F00") |
Diese Version von dqrng ist noch nicht auf CRAN verfügbar, kann aber über drat installiert werden:
|
1 2 3 4 |
if (!requireNamespace("drat", quietly = TRUE)) install.packages("drat") drat::addRepo("daqana") install.packages("dqrng") |
Insgesamt werden wir sehen, dass es möglich ist schneller als R Zufallsstichproben zu erstellen ohne Kompromisse bei der statistischen Qualität eingehen zu müssen. Feedback ist immer willkommen!
Benchmarks
Die folgenden Benchmarks wurden mit R Version 3.5.3 (2019-03-11) unter Debian GNU/Linux 9 (stretch) auf einem Intel® Core™ i7-6600U CPU @ 2.60GHz erstellt. YMMV, wie immer bei Benchmarks …
Durch die Kombination schneller Zufallszahlgeneratoren mit schnellen Methoden zur Erzeugung ganzer Zahlen erhält man gute Performance beim Ziehen mit Zurücklegen:
|
1 2 3 4 5 6 7 8 9 10 |
library(dqrng) m <- 1e6 n <- 1e4 bm <- bench::mark(sample.int(m, n, replace = TRUE), sample.int(1e4*m, n, replace = TRUE), dqsample.int(m, n, replace = TRUE), dqsample.int(1e4*m, n, replace = TRUE), check = FALSE) knitr::kable(bm[, 1:6]) |
| expression | min | mean | median | max | itr/sec |
|---|---|---|---|---|---|
| sample.int(m, n, replace = TRUE) | 94.9µs | 112.5µs | 112.2µs | 3.78ms | 8886.863 |
| sample.int(10000 * m, n, replace = TRUE) | 170.9µs | 202.8µs | 200.9µs | 3.08ms | 4930.094 |
| dqsample.int(m, n, replace = TRUE) | 31.5µs | 36.2µs | 35.7µs | 183.63µs | 27596.508 |
| dqsample.int(10000 * m, n, replace = TRUE) | 35.7µs | 41.3µs | 39.8µs | 3.06ms | 24186.338 |
|
1 2 |
plot(bm) |
Beachte dass für 10^10 Zahlen der “long-vector support” in R benötigt wird.
Beim Ziehen ohne Zurücklegen muss man sich Gedanken machen wie man sicher stellt, dass kein Eintrag sich wiederholt. Wenn mehr als 50% der Gesamtpopulation gezogen werden sollen, so mischt dqrng die Gesamtpopulation bis zur benötigten Größe. R nutzen eine sehr ähnlichen Algorithmus, jedoch ist dqrng auch hier schneller:
|
1 2 3 4 5 6 7 8 |
library(dqrng) m <- 1e6 n <- 6e5 bm <- bench::mark(sample.int(m, n), dqsample.int(m, n), check = FALSE, min_iterations = 50) knitr::kable(bm[, 1:6]) |
| expression | min | mean | median | max | itr/sec |
|---|---|---|---|---|---|
| sample.int(m, n) | 11.91ms | 14.09ms | 13.02ms | 26.11ms | 70.98544 |
| dqsample.int(m, n) | 6.93ms | 8.22ms | 8.17ms | 9.99ms | 121.65317 |
|
1 2 |
plot(bm) |
Wird ein geringerer Teil der Population gezogen, so verwendet dqrng die Verwerfungsmethode mit einem Set als Speicher der bereits gezogenen Element. Grundsätzlich kann R einen ganz ähnlichen Algorithmus auf Basis eines Hashsets verwenden, tut dies standardmäßig aber erst für größere Populationen, obwohl es schon früher schneller als die Standardmethode ist. Der in dqrng genutzte und auf einem Bitset beruhende Algorithmus ist dabei noch schneller:
|
1 2 3 4 5 6 7 8 9 |
library(dqrng) m <- 1e6 n <- 1e4 bm <- bench::mark(sample.int(m, n), sample.int(m, n, useHash = TRUE), dqsample.int(m, n), check = FALSE) knitr::kable(bm[, 1:6]) |
| expression | min | mean | median | max | itr/sec |
|---|---|---|---|---|---|
| sample.int(m, n) | 788µs | 1.72ms | 1.87ms | 5.01ms | 583.0849 |
| sample.int(m, n, useHash = TRUE) | 229µs | 269.6µs | 266.15µs | 453.71µs | 3709.2196 |
| dqsample.int(m, n) | 113µs | 130.4µs | 129.49µs | 2.98ms | 7668.6069 |
|
1 2 |
plot(bm) |
Wird der Stichprobenanteil weiter reduziert, so verwendet auch dqrng ein Hashset. Die beiden Hashset basierten Algorithmen sind vergleichbar schnell und deutlich schneller als die Standardmethode in R:
|
1 2 3 4 5 6 7 8 9 |
library(dqrng) m <- 1e6 n <- 1e2 bm <- bench::mark(sample.int(m, n), sample.int(m, n, useHash = TRUE), dqsample.int(m, n), check = FALSE) knitr::kable(bm[, 1:6]) |
| expression | min | mean | median | max | itr/sec |
|---|---|---|---|---|---|
| sample.int(m, n) | 452.72µs | 1.29ms | 1.56ms | 4.69ms | 777.0096 |
| sample.int(m, n, useHash = TRUE) | 3.98µs | 5.38µs | 5.03µs | 62.94µs | 185978.0662 |
| dqsample.int(m, n) | 3.59µs | 4.14µs | 4.05µs | 26.85µs | 241454.9528 |
|
1 2 |
plot(bm) |
Für größere Gesamtpopulationen verwendet R standardmäßig ein Hashset. Jedoch ist die Version in dqrng schneller:
|
1 2 3 4 5 6 7 8 |
library(dqrng) m <- 1e10 n <- 1e5 bm <- bench::mark(sample.int(m, n), dqsample.int(m, n), check = FALSE) knitr::kable(bm[, 1:6]) |
| expression | min | mean | median | max | itr/sec |
|---|---|---|---|---|---|
| sample.int(m, n) | 5.3ms | 5.79ms | 5.48ms | 9.5ms | 172.7915 |
| dqsample.int(m, n) | 1.7ms | 2.04ms | 1.97ms | 5.08ms | 490.8562 |
|
1 2 |
plot(bm) |
Details
Die folgenden Methoden werden beim Ziehen ohne Zurücklegen verwendet. Die Algorithmen werden in R-artigem Pseudocode dargestellt, obwohl die tatsächliche Umsetzung in C++ erfolgte. Für Stichprobenanteile über 50% wird partiell nach Fisher-Yates gemischt:
|
1 2 3 4 5 6 7 |
no_replace_shuffle <- function(m, n) { tmp <- seq_len(m) for (i in seq_len(n)) swap(tmp[i], tmp[i + random_int(m-i)]) tmp[1:n] } |
wobei random_int(m-i) eine Zufallszahl aus [0, m-i] liefert. Nachdem hier die Gesamtpopulation im Speicher gehalten werden muss, ist diese Methode nur für hohe Stichprobenanteile geeignet. Man könnte erwarten, dass Reservoir Sampling bei geringeren Anteilen gut funktioniert. In meinen Tests waren jedoch Set-basierte Algorithmen schneller:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
no_replace_set <- function(m, n) { result <- vector(mode = "...", length = n) # integer or numeric elems <- new(set, m, n) # set object for storing n objects out of m possible values for (i in seq_len(n)) while (TRUE) { v = random_int(m) if (elems.insert(v)) { result[i] = v break } } result } |
Dabei gibt elems.insert(v) den Wert TRUE zurück, wenn v nicht Teil von elems war und daher eingefügt wurde, andernfalls wird FALSE zurückgegeben. Es gibt verschiedene Möglichkeiten solch ein Set umzusetzen. Für mittlere Stichprobenanteile (derzeit zwischen 0.1% und 50%) verwendet dqrng ein Bitset, d.h. einen Vektor aus m Bits von denen jeder einen der möglichen Werte repräsentiert. Für niedrigere Stichprobenanteile wird der Speicherbedarf des Bitsets zum Problem, weshalb dann ein Hashset zum Einsatz kommt.1 Hier skaliert der Speicherbedarf nur mit n und nicht mit m. Man könnte erwarten, dass der auf Robert Floyd zurückgehende Algorithmus besser sein sollte, aber das war in meinen Tests nicht der Fall. Vermutlich weil das Ergebnis nochmals gemischt werden muss um eine zufällige _Permutation_ und nicht nur eine zufällige _Kombination_ zu erhalten.
- Für die Spezialisten: Open addressing mit power-of-two Größe zwischen 1.5 und 3 mal
n, identity hash und quadratic probing.