Der häufigste Fehler bei der Anwendung der Funktion sample()

Die Funktion sample() wird verwendet, um Stichproben zu erzeugen. Damit dies nicht zu unerwünschtem Verhalten führt, muss man wissen, dass der Aufruf an sample.int() weitergereicht wird, wenn die Menge, aus der ein Objekt ausgewählt werden soll, nur ein Element besitzt. Es wird ein Beispiel ausführlich besprochen, bei dem eine naheliegende Implementierung zu unerwarteten Ergebnissen führt.

Einordnung des Artikels

Einführung

Die Funktion sample() aus dem Paket base erlaubt es das Ziehen einer Stichprobe zu simulieren. Dazu kann man vorgeben:

sample(x, size, replace = FALSE, prob = NULL)

Da es für die Funktion sample() eine "interger"-Version gibt, nämlich sample.int(), die die Grundmenge x anders konfiguriert, kann es leicht zu einem Fehler kommen, wenn zwar im Quelltext ein Aufruf von sample() steht, intern aber sample.int() aufgerufen wird.

sample.int(n, size = n, replace = FALSE, prob = NULL,
useHash = (!replace && is.null(prob) && size <= n/2 && n > 1e7))

Sämtliche Eingabewerte von sample.int() werden hier nicht besprochen, die Fehlerquelle liegt im Argument n. Dazu wird im Folgenden ein Beispiel (spezieller Zufallsweg auf den Kanten eines Würfels) ausführlich behandelt, das den Fehler aufzeigt. Das Beispiel ist so gewählt, dass der dabei auftretende "Fehler" unverständlich ist, wenn man nur die Funktion sample() und nicht sample.int() kennt. Anschließend wird gezeigt, wie man den Fehler vermeiden kann.

Beispiel: Ein spezieller Zufallsweg auf den Kanten eines Würfels

Die Aufgabenstellung

Abbildung 1 zeigt oben den Einheitswürfel, dessen Ecken mit den Koordinaten (0, 0, 0), ..., (1, 1, 1) beschriftet sind. Es soll ein Zufallsweg von (0, 0, 0) nach (1, 1, 1) bestimmt werden, der folgende Eigenschaften besitzt:

In Abbildung 1 unten ist der Weg

(0, 0, 0) → (0, 1, 0) → (0, 1, 1) → (1, 1, 1)

eingezeichnet, der diese Bedingungen erfüllt.

Abbildung 1: Oben: Auf dem dreidimensionalen Einheits-Würfel sollen Zufallswege bestimmt werden, die von (0, 0, 0) nach (1, 1, 1) führen. Dabei muss der Weg entlang Kanten laufen und er darf nur aus drei Kanten zusammengesetzt sein. Unten: von den 6 möglichen Wegen ist ein spezieller Weg dargestellt (orange).Abbildung 1: Oben: Auf dem dreidimensionalen Einheits-Würfel sollen Zufallswege bestimmt werden, die von (0, 0, 0) nach (1, 1, 1) führen. Dabei muss der Weg entlang Kanten laufen und er darf nur aus drei Kanten zusammengesetzt sein. Unten: von den 6 möglichen Wegen ist ein spezieller Weg dargestellt (orange).

Man kann die Aufgabenstellung leicht verallgemeinern und den Weg auf einem n-dimensionalen Würfel suchen; die folgenden Quelltexte sind so gestaltet, dass man diese Verallgemeinerung leicht vornehmen kann. Da hier lediglich das Verhalten der Funktion sample() besprochen werden soll, wird nur der dreidimensionale Würfel betrachtet.

Eine erste (fehlerfreie) Implementierung mit sample()

Da man einen kürzesten Weg wählen muss, besteht er aus drei aufeinanderfolgenden Schritten, wobei man sich in jedem Schritt in einer anderen Richtung bewegen muss. Kürzt man die x-, y, z-Richtung durch die Zahlen 1, 2, 3 ab, so wird ein Weg von (0, 0, 0) nach (1, 1, 1) durch die Angabe der drei Zahlen festgelegt. Der Weg

(0, 0, 0) → (0, 1, 0) → (0, 1, 1) → (1, 1, 1)

kann dann durch die Abfolge

2, 3, 1

symbolisiert werden. Damit ist auch klar, dass es insgesamt 3 · 2 · 1 = 6 mögliche Wege gibt. (In n Dimensionen gibt es n! mögliche Wege.)

Und es ist sofort klar, durch welches Zufallsexperiment die Auswahl eines Weges realisiert wird: Man zieht aus den Zahlen (1, 2, 3) nacheinander drei Zahlen ohne Zurücklegen.

Die Implementierung besteht dann nur noch aus einem Aufruf der Funktion sample(), die geeignet konfiguriert wird:

set.seed(17)
sample(x = (1:3), size = 3)
# 2 3 1

Für x wird die Menge der Zahlen von 1 bis 3 angegeben und es werden 3 Zahlen gezogen. Weitere Argumente muss man nicht setzen, da hier die default-Werte das gewünschte Verhalten liefern.

Der Nachteil dieser Implementierung ist offensichtlich: Erzeugt werden Richtungen; ist man tatsächlich an den Koordinaten der Punkte interessiert, die der Weg passiert, muss man die Richtungen erst umrechnen

Eine zweite (fehlerhafte) Implementierung mit sample()

Die folgende Implementierung versucht die Abfolge der Koordinaten zu erzeugen, die entlang des Zufallsweges passiert werden. Dazu startet man im Ursprung, der durch den Vektor c(0, 0, 0) modelliert wird (Zeile 1). Der Zufallsweg ist vollständig, wenn man den Punkt c(1, 1, 1) erreicht; entsprechend ist die while-Schleife konfiguriert (Zeile 3 bis 8).

In jedem Iterationsschritt wird aus den Koordinaten, die gleich 0 sind, eine Koordinate ausgewählt (Zeile 4) und gleich 1 gesetzt (Zeile 5). Man könnte die beiden Anweisungen in Zeile 4 und 5 zu einer Anweisung zusammenfassen, dies geschieht hier nicht, um den Fehler der Implementierung besser zeigen zu können.

Nach 3 Iterationsschritten landet man im Punkt c(1, 1, 1) und hat dabei die Koordinaten der Zwischenstationen berechnet (für eine hochdimensionale Anwendung mag diese Implementierung günstiger erscheinen als die erste.)

v <- c(0, 0, 0)
print(v)
while(!identical(v, c(1, 1, 1))){
  idx <- which(v == 0)
  idx.0 <- sample(x = idx, size = 1)
  v[idx.0] <- 1
  print(v)
}

Führt man das Skript aus mit set.seed(17) , erhält man die Ausgaben:

# [1] 0 0 0
# [1] 0 1 0
# [1] 0 1 1
# [1] 1 1 1

Dagegen erhält man mit set.seed(1) die Folge von Punkten:

# [1] 0 0 0
# [1] 1 0 0
# [1] 1 0 1
# [1] 1 0 1
# [1] 1 0 1
# [1] 1 1 1

Wie kann es sein, dass im zweiten Test aus dem Punkt 1 0 1 (zum ersten Mal in Zeile 3) die 0 erst im dritten Anlauf entfernt wird?

Dazu wird nachvollzogen, welche Berechnungen bei set.seed(1) (zweiter Test) ausgeführt werden (Abbildung 2 versucht den Ablauf der Berechnungen darzustellen):

Die letzte Erwartung ist falsch und das Verhalten von sample() wird nur verständlich, wenn man weiß, dass der Aufruf an sample.int() weitergereicht wird. In diesem Fall (mit idx gleich 2) wird somit aufgerufen:

sample.int(n = 2, size = 1)

Abbildung 2: Darstellung des Durchlaufs der while-Schleife aus dem Skript zur Berechnung eines Zufallsweges von (0, 0, 0) nach (1, 1, 1). Sobald idx ein Objekt der Länge 1 ist, wird der Aufruf von sample() an sample.int() weitergereicht, was hier zu einem unerwünschten Ergebnis führt.Abbildung 2: Darstellung des Durchlaufs der while-Schleife aus dem Skript zur Berechnung eines Zufallsweges von (0, 0, 0) nach (1, 1, 1). Sobald idx ein Objekt der Länge 1 ist, wird der Aufruf von sample() an sample.int() weitergereicht, was hier zu einem unerwünschten Ergebnis führt.

Die Implementierung von sample.int() liest den Eingabewert n = 2 als die Abkürzung für (1:2) , also als zweielementige Menge. Und dieser Aufruf erzeugt entweder 1 oder 2 (jeweils mit Wahrscheinlichkeit 1/2). Wird aber 1 erzeugt, so wird die erste Komponente von v nochmals gleich 1 gesetzt (und es entsteht Zeile 4 der Ausgabe). Erst wenn die 2 erzeugt wird, erhält man die Ausgabe aus Zeile 6.

Die Fehlerquelle

In der zweiten Implementierung des Zufallsexperimentes wurde nicht berücksichtigt:

Hat das Objekt x im Aufruf von sample() die Länge 1, so wird der Aufruf an sample.int() weitergereicht.

Es wird dann aufgerufen:

sample.int(n = x, size = 1)

und somit ein Element aus der Menge (1:x) ausgewählt. Dass das erste Argument beim Aufruf von sample() tatsächlich mit x = idx benannt wurde und nicht nur idx.0 <- sample(idx, size = 1) geschrieben wurde, verhindert dieses Weiterreichen nicht.

Umgekehrt sollte man bei der Anwendung der Funktion sample() immer bedenken: Wird an sample() eine Menge x übergeben, so kann der Aufruf unerwünschtes Verhalten erzeugen, wenn die Menge nur aus einem Element besteht. Das obige Beispiel wurde gerade so konstruiert, um diese Fehlerquelle aufzuzeigen (und auch, dass der Fehler nicht in jedem Durchlauf auftreten muss und somit beim Testen leicht unentdeckt bleiben kann).

Eine sichere Implementierung von sample()

Das folgende Skript zeigt eine Funktion safeSample(), die das oben beschriebene unerwünschte Verhalten von sample() vermeidet. Es wird nicht behauptet, dass man diese sichere Implementierung immer einsetzen sollte, wenn man mit sample() arbeiten möchte. Aber man sollte diese (oder eine ähnliche) Implementierung verwenden, wenn der oben beschriebene Fehler auftreten kann.

safeSample <- function(x, size, replace = FALSE, prob = NULL){
  if(identical(length(x), 1L)){
    return(x)
  } else {
    return(sample(x = x, size = size, replace = replace, prob = prob))
  }
}

Wenn die Menge x nur ein Element hat, muss dieses Element zurückgegeben werden (Zeile 2 und 3). Andernfalls wird der Aufruf an sample() weitergereicht.

Man kann die Implementierung auch erweitern und sich für den Fall identical(length(x), 0L) eine geeignete Rückgabe (oder Fehlermeldung) überlegen.

Mit Hilfe der Funktion safeSample() kann jetzt die zweite Implementierung verbessert werden: an die Stelle des Aufrufs von sample() tritt der Aufruf von safeSample():

set.seed(1)
v <- c(0, 0, 0)
print(v)
while(!identical(v, c(1, 1, 1))){
  idx <- which(v == 0)
  idx.0 <- safeSample(x = idx, size = 1)
  v[idx.0] <- 1
  print(v)
}
# [1] 0 0 0
# [1] 1 0 0
# [1] 1 0 1
# [1] 1 1 1

Mit set.seed(17) erhält man die Ausgabe, die oben schon gezeigt wurde.