Archiv des Autors: Ralf Stubner

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.

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.

dqsample: Eine faire Alternative zu ‚base::sample‘

Das Erzeugen von zufälligen Stichproben eines Datensatzes wird in der Statistik oder Data Science oft angewendet. Für diese Aufgabe gibt es in R die Funktion base::sample. Leider ist der dabei verwendete Algorithmus nicht vollständig fair. Dies wurde zuletzt auf R-devel diskutiert, was auch die Motivation für das dqsample Paket darstellt. Derzeit ist dqsample nicht auf CRAN, kann aber über drat installiert werden:

Beispiel

Wählt man viele Zahlen zufällig aus, so sollten die Dichten der geraden und ungerade Zahlen in etwa gleich und konstant sein. Mit base::sample ist dies jedoch nicht er Fall:

plot of chunk base

Oder mit verändertem Parameter:

plot of chunk base-oszi

Dieses spezielle Beispiel stammt von Duncan Murdoch.

Daniel Lemire (2018, <arXiv:1805.1094>) hat einen alternativen Algorithmus vorgeschlagen, der in dqsample verwendet wird. Mit diesem Algorithmus gibt es keine Asymmetrie zwischen geraden und ungeraden Zahlen:

plot of chunk dqsample

Ursachenforschung

Intern benötigt die base::sample() Funktion zufällige ganze Zahlen die gleichmäßig in dem halb-offenen Bereich [0, n) verteilt sind.Um diese zu erzeugen, verwendet R zufällige Gleitkommazahlen aus [0, 1), multipliziert mit n und schneidet den Nachkommateil ab. Mit reellen Zahlen im mathematischen Sinn wäre das korrekt. Aber das ist hier nicht der Fall.

Standardmäßig verwendet R einen 32 bit Mersenne-Twister zum Erzeugen von Zufallszahlen. Dieser erzeugt ganze Zahlen aus [0, 2^32), die durch 2^32geteilt werden um Gleitkommazahlen aus [0, 1) zu erhalten.Wenn wir das oben beschriebene Verfahren invertieren, so können wir sehen wie viele mögliche Zufallszahlen zu einem Ergebnis gehören.  Mit sample(6, 10, replace = TRUE) kann man zum Beispiel den Wurf von zehn Würfeln simulieren. Da 2^32 kein Vielfaches von sechs ist, kann die Verteilung aber nicht gleichverteilt sein:

Eins und vier haben eine minimal geringere Wahrscheinlichkeit als die übrigen Zahlen. Dieser Effekt steigt mit der Größe des Datensatzes, aus dem gewählt wird. Verwenden wir das oben definierte m so können wir sehen, woher die ungleiche Verteilung gerader und ungerader Zahlen herkommt:

Während jeweils zwei Zahlen auf eine der ungeraden Zahlen abgebildet werden, sind es drei Zahlen für jede gerade Zahl. Dieses Muster verschiebt sich bei der Hälfte der möglichen Ergebnisse, was das erste Bild oben erklärt. Das Muster verschiebt sich öfter, wenn man sich von m entfernt, was die Oszillationen im zweiten Bild erklärt. Irgendwann werden die Oszillationen zu schnell, um sie in der Dichte sehen zu können. Das bedeutet aber nicht, dass die Verteilung fair ist. Für m - 2^20verschiebt sich das Muster beispielsweise zwischen 982 und 983:

Vorher sind gerade Zahlen wahrscheinlicher als ungerade Zahlen. Danach ist es umgekehrt.

Zusammenfassung

Der von base::sample verwendet Algorithmus ist nicht fair, was zusammen mit großen Datenmengen einen merkbaren Effekt haben kann. Das dqsample Paket stellt für die wichtigsten Anwendungsfälle einen fairen Algorithmus bereit. Es kann als direkter Ersatz für base::sample verwendet werden.

tikzDevice hat ein neues zu Hause

Im Februar diesen Jahres wurde das tikzDevice Paket auf CRAN auf ORPHANED gestellt. In der Folge suchten Kirill Müller und Yihui Xie eine neuen Maintainer. Als wir das einige Zeit später lasen, beschlossen wir uns hier einzubringen. Nach einem kurzen E-Mail-Dialog mit Yihui Xie wurde das GitHub Repository auf unsere Organization übertragen: https://github.com/daqana/tikzDevice. Derweil haben wir die bestehenden Warnungen aus den CRAN checks behoben und Version 0.12 hochgeladen.

Was kann man mit tikzDevice tun? es is ein graphics device vergleichbar zu pdf() oder png(). Aber an Stelle einer Bilddatei, die in einen Report als externe Graphik eingefügt werden kann, erzeugt es Dateien im TikZ format durch die LaTeX die Graphiken erstellt. Das ermöglicht konsistente Schriften für Text und Graphiken sowie die Verwendung von TeXs Fähigkeiten zur Darstellung mathematischer Formeln. Die PDF-Vignette enthält viele Beispiele.

Man kann es sogar in einem R-Markdown Dokument nutzen. Ein Dokument mit

im YAML Kopf und

in einem setup Block wird für Text und Graphiken die gleiche Schrift nutzen, wenn letztere mit dev = "tikz" erzeugt werden. In diesem Beispiel Palatino mit Minuskelziffern:

Palatino mit Minuskelziffern