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

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Wahrscheinlichkeit und Statistik
      • Wahrscheinlichkeitsverteilungen in R
      • Die Funktion sample() zum Erzeugen von Stichproben zu selbstdefinierten diskreten Zufallsvariablen
      • Der häufigste Fehler bei der Anwendung der Funktion sample()

Einführung

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

• Aus welcher Menge x die Stichprobenelemente gewählt werden,
• wie groß die Stichprobe ist (size),
• mit welchen Wahrscheinlichkeiten die Elemente aus x gezogen werden (prob) und
• ob das Ziehen mit Zurücklegen (replace = TRUE ) oder ohne Zurücklegen (replace = FALSE ) realisiert werden soll.

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:

• man darf sich nur auf den Kanten bewegen und
• es muss ein kürzester Weg gewählt werden.

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).

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):

• Im ersten Durchlauf der while-Schleife besteht idx aus den drei Zahlen (1, 2, 3). Daraus wird mit sample() eine Zahl ausgewählt, hier die 1, wodurch die erste Komponente von v gleich 1 gesetzt wird (Ausgabe in Zeile 2).
• Im zweiten Durchlauf der while-Schleife besteht idx aus den beiden Zahlen (2, 3); jetzt wird die 3 ausgewählt und es entsteht der Vektor v = (1, 0, 1), siehe Zeile 3 der Ausgabe.
• Im dritten Durchlauf der while-Schleife besteht idx nur noch aus der Zahl 2 (und soll – so ist die Implementierung beabsichtigt – als einelementige Menge interpretiert werden). Aber hier geschieht offensichtlich etwas Unerwartetes: Man erwartet, dass durch sample() die Zahl 2 gewählt werden muss und somit sofort v = (1, 1, 1) entsteht.

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.

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.