Schnellere Zufallsstichproben mit dqrng

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:

Diese Version von dqrng ist noch nicht auf CRAN verfügbar, kann aber über drat installiert werden:

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:

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

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:

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

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:

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

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:

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

Für größere Gesamtpopulationen verwendet R standardmäßig ein Hashset. Jedoch ist die Version in dqrng schneller:

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

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:

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:

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.

  1. Für die Spezialisten: Open addressing mit power-of-two Größe zwischen 1.5 und 3 mal n, identity hash und quadratic probing.

Eine kleine Essensumfrage – Ergebnisse

Ende letzten Jahres haben wir eine kleine (nicht-repräsentative) Essensumfrage gestartet: Welche Essen sind besonders beliebt, gibt es verschiedene “Esstypen” und korrelieren Essvorlieben mit Reisezielen?

Die Ergebnisse

Nun präsentieren wir die Ergebnisse. Es haben 60 Personen teilgenommen, von denen die meisten 20 – 40 Jahre alt waren. Ein erstaunlich hoher Anteil von über 15 % gab an, sich vegan zu ernähren. Das Verhältnis zwischen weiblichen und männlichen Teilnehmenden war ausgeglichen. Im Mittel wurden etwa 60% der abgefragten Lebensmittels gerne gegessen.

Körnerbrötchen, Orangen und Schokolade gehörten zu den beliebtesten Essensgegenständen. Alge hatten immerhin 17 % noch nie gegessen. Bei einigen der abgefragten Essensgegenständen hat sich eine Abneigung im Laufe des Lebens in eine Zuneigung gewandelt. Hier wurden auffällig viele Gemüse genannt, allen voran Oliven. Insgesamt mochten vegan lebende Menschen mehr der abgefragen Essen als die restlichen Teilnehmenden.

Schließlich gab es gewisse Zusammenhänge zwischen den Essensvorlieben und demographischen Merkmalen, wie besuchten Kontinenten, Lieblingsfarben und genutzten Betriebssystemen. Beispielsweise zeigte sich, dass Inwger besonders bei den Teilnehmenden beliebt war, die auch Nordamerika besucht hatten. Andererseits waren z.B. weiße Brötchen besonders bei Südamerika-Reisenden unbeliebt. Weitere Ergebnisse und Erklärungen gibt es über die Links.

Bis bald und liebe Grüße!

dqrng v0.1.0: breaking changes

Eine neue Version des dqrng Pakets ist auf den CRAN Servern verfügbar. Diese Version bringt zwei Änderungen des bisherigen Verhaltens:

  • Ein int-Vektor wird an Stelle eines einzelnen int als Seed verwendet (Aaron Lun in #10). Verwendung eines einzelnen int führt zu einem anderen Zustand als zuvor.

  • Der Mersenne-Twister Generator wurde entfernt. Xoroshiro128+ ist nun der default.

Die erste Änderung ist durch den Wunsch motiviert mehr als 32 bits an Zufall
als Seed verwenden zu können. Damit wurde das zuvor genutzte “scrambling” des einzelnen int überflüssig. Die neue Methode generateSeedVectors() nutzt R's RNG um eine Liste an int-Vektoren zu erzeugen, die als Seed verwendet werden könnnen.

Die zweite Änderung bezieht sich auf einen Satz aus der Anleitung zu R: Nor should the C++11 random number library be used …. I vermute, dass sich das auf die impementierungs-abhängigen Verteilungen und nicht die Generatoren bezieht, aber grundsätzloch sollte man WRE möglichst genau folgen. Also muss std::mt19937_64 gehen und kann wegen eines nicht gemerged pull requests nicht durch boost::random::mt19937_64 ersetzt werden. Statt eine gefixte Versio auszuliefern entschied ich mich MT zu entfernen da:

  • MT einige statistische Tests nicht besteht.
  • MT langsammer ist als die anderen Generatoren.
  • MT der einzige Generator war, der keine Unterstützung für parallele Verarbeitung mitbringt.

Daneben wurde noch std::random_device aus random von C++11 während der Initialisierung genutzt. Dies ist nun unnötig, da der initiale Zustand des default RNG über die Techniken von generateSeedVectors() bestimmt wird.

dqrng v0.0.5: Neue und aktualisierte RNGs

Eine neue Version des dqrng Pakets ist auf den CRAN Servern verfügbar. Vielen Dank an die CRAN Maintainer im Allgemeinen und Uwe Ligges im Besonderen für ihren Einsatz.

In dieser Version gibt es einen neuen Zufallszahlengenerator, der zusammen mit den Verteilungsfunktionen benutzt werden kann: Die 64 bit Version des Threefry Engines mit 20 Runden (Salmon et al., 2011 <doi:10.1145/2063384.2063405>, http://www.thesalmons.org/john/random123/) aus dem sitmo Paket. Vielen Dank an James Balamuta für das Bereitstellen der C++11 Version von Threefry in seinem Paket.

Zusätzlich wurden die PCG-Header auf die neueste Version aktualisiert. Vielen Dank an Aaron Lun für die Erinnerung und den Hinweis auf zwei BioConductor-Pakete die dgrng nutzen: DropletUtils und scran.

Daneben gab es noch zwei technische Änderungen: Einerseits ist das Paket nun vorbereitet für die Zukunft, nachdem -DSTRICT_R_HEADERS in src/Makevars zu PKG_CPPFLAGS hinzugefügt wird.

Andererseits gibt es nun Unit-Tests für die C++ Schnittstelle. Zunächst hatte ich Catch probiert, nachdem ich bereits testthat nutze. Allerdings gefiel es mir nicht Testcode in die endgültige Bibliothek zu integrieren. Stattdessen gibt es nun C++ Dateien, die während des Tests durch Rcpp::sourceCpp() kompiliert und danach ausgeführt werden. Das ist zudem ein realistischerer Test der C++ Schnittstelle.

Eine kleine Essensumfrage

In diesem Beitrag geht es um ein kleines Seitenprojekt: Eine (nicht-repräsentative) Essensumfrage um einen Eindruck zu bekommen, welches Essen wie (un-)beliebt ist, ob Essvorlieben mit Reisezielen korrelieren, ob es verschiedene “Esstypen” gibt und vieles mehr.

Die Idee

Während des Mittags reden wir oft über Essen. Über verschiedene Geschmäcker, wie (un-)beliebt bestimmte Essen sind und ob es vielleicht verschiedene Esstypen gibt (Leute, die Lakritze mögen, essen ungern XY?). Um etwas Licht in diese Fragen zu bringen, weil wir keine zugänglichen Studien finden konnten, und mit allen Werkzeugen der Web-Entwicklung und Statistik im Handgepäck, entschlossen wir uns für ein Seitenprojekt.

Nun präsentieren wir eine kleine Essensumfrage! Da wir keine repräsentative Stichprobe auswählen, lassen sich die Ergebnisse nicht verallgemeinern. Die Ergebnisse sagen nur etwas über die teilnehmenden Personen aus, nicht über Deutschland, sicher nicht über die ganze Welt. Wir hoffen, dass die Ergebnisse uns und dir trotzdem Freude bereiten und einige interessante Zusammenhänge liefern.

Weitere Infos

Nachdem du dich durch drei mittelkurze Seiten mit Fragen geklickt hast, gibt es außerdem einige vorläufige Ergebnisse zu sehen und wie die Anderen im Vergleich zu dir geantwortet haben. Außerdem gilt selbstverständlich:

  • Alle Ergebnisse sind vollständig anonym und werden nicht kommerziell genutzt.
  • Über den Verlauf unserer kleinen Umfrage werden wir hier in diesem Blog-Beitrag informieren.
  • Wenn du die Umfrage weitergeben möchtest, bist du herzlich dazu eingeladen.

Wir freuen uns über Kommentare und Fragen. Viel Spaß!