Die Funktion sample() zum Erzeugen von Stichproben zu selbstdefinierten diskreten Zufallsvariablen

Die Funktion sample() wird verwendet, um Stichproben zu erzeugen. Sie lässt sich so konfigurieren, dass man die Wahrscheinlichkeitsverteilungen von beliebigen selbstdefinierten diskreten Zufallsvariablen einsetzen kann. Zudem kann man das Ziehen mit beziehungsweise ohne Zurücklegen realisieren.

walter

02 Dez. 2020

Wahrscheinlichkeitsverteilung Ziehen ohne Zurücklegen Fixpunkt diskrete Zufallsvariable Permutation replicate() Zufallsvariable Stichprobe Ziehen mit Zurücklegen Dataframe sample()

Inhaltsverzeichnis

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

Zum besseren Verständnis werden Kenntnisse über Zufallsvariablen vorausgesetzt, wie sie etwa in Grundbegriffe der Wahrscheinlichkeitsrechnung: Diskrete und stetige Zufallsvariablen vermittelt werden.

Ein typischer Fehler in der Anwendung der Funktion sample() wird in Der häufigste Fehler bei der Anwendung der Funktion sample() besprochen. Er beruht darauf, dass der Aufruf von sample() an die Funktion sample.int() weitergereicht wird, wenn die Stichprobe aus einer einelementigen Menge gezogen wird – und dann zu unerwartetem Verhalten von sample() führt. Auf den ersten Blick ist es kaum vorstellbar, wie diese Situation eintreten kann, dort wird ein Beispiel besprochen, das zeigt, dass man sehr leicht in diese Falle tappen kann.

Einführung

Die Basis-Pakete von R erlauben den Umgang mit einer Vielzahl von Wahrscheinlichkeitsverteilungen; einige davon wurden in Wahrscheinlichkeitsverteilungen in R besprochen, weitere findet man unter

?Distributions

Zumindest für Verteilungen von Zufallsvariablen X, die endlich viele Werte annehmen, wird eine einfache Möglichkeit vorgestellt, wie man sie mit Hilfe eines Dataframes modellieren kann. Auf Zufallsvariablen mit abzählbar vielen Werten (wie etwa die Poisson-Verteilung) kann man dieses Modell nicht übertragen.

Möchte man eine Simulation eines Zufallsexperimentes mit dieser Zufallsvariable X durchführen – also Zufallszahlen erzeugen, die gemäß der Verteilung von X "ausgewürfelt" werden –, kann man die Funktion sample() einsetzen. So kann man auch mit Verteilungen arbeiten, die nicht in ?Distributions enthalten sind.

Es soll ausdrücklich betont werden, dass die Modellierung einer Zufallsvariable als Dataframe nicht zwingend erforderlich ist, um die Funktion sample() anzuwenden. Sie ist aber dann zu empfehlen, wenn die Zufallsvariable mehrmals benötigt wird und nicht nur eine Stichprobe erzeugt werden soll.

In diesem Kapitel wird:

zunächst die oben erwähnte Modellierung einer Zufallsvariablen durch ein Dataframe beschrieben,
dann die Funktion sample() und vor allem die Bedeutung ihrer Eingabewerte vorgestellt.
Dazu werden typische Anwendungsbeispiele gezeigt.
Insbesondere wird gezeigt, wie man Simulationen mit Hilfe von replicate() wiederholen kann, ohne Schleifen zu implementieren.

Modellierung einer Zufallsvariable als Dataframe

Die hier beschriebene Modellierung einer Zufallsvariable durch ein Dataframe ist identisch mit der Beschreibung in Grundbegriffe der Wahrscheinlichkeitsrechnung: Die Zufallsvariable. Dort wird auch der mathematische Hintergrund erläutert, insbesondere der Zusammenhang zwischen Zufallsvariablen und einer Wahrscheinlichkeitsverteilungen.

Um mit Zufallsvariablen und ihren Verteilungen zu arbeiten, benötigt man ein möglichst einfaches R-Objekt, in dem die Eigenschaften der Zufallsvariable abgelegt werden. Man kann dann – je nach gewünschter Anwendung – leicht darauf zugreifen und geeignete Service-Funktionen selber implementieren.

Es wird hier ein einfacher Zugang gewählt, der immer dann anwendbar ist, wenn eine Zufallsvariable eine Wertemenge mit endlich vielen reellen Werten besitzt: Die Zufallsvariable wird durch ein Dataframe modelliert, das zwei Spalten besitzt:

Eine Spalte value für die möglichen Werte der Zufallsvariable.
Eine Spalte prob für die Wahrscheinlichkeiten der entsprechenden Werte.

Es ist klar, dass die Summe aller Wahrscheinlichkeiten eins ergeben muss, dies und weitere Anforderungen an das Dataframe werden weiter unten diskutiert.

Die Eigenschaften von Dataframes wurden in Dataframes in R: der Datentyp data frame und Dataframes in R: Anwendungen ausführlich erklärt; die nötigen Kenntnisse darüber werden hier vorausgesetzt.

Erzeugen von Zufallsvariablen

Das folgende Skript zeigt, wie eine Zufallsvariable als Dataframe erzeugt wird, die eine gezinkte Münze beschreibt. Dabei werden die Ergebnisse des Zufallsexperimentes durch 0 und 1 modelliert (anstelle von Kopf und Zahl); die Wahrscheinlichkeiten dafür, dass die Münze 0 oder 1 zeigt, sollen 0.4 beziehungsweise 0.6 betragen:

values <- c(0, 1)
probs <- c(0.4, 0.6)

coin <- data.frame(value = values, prob = probs)
print(coin, digits = 2)
#   value  prob
# 1     0 0.40
# 2     2 0.60

Zeile 1 und 2: Die Werte und die entsprechenden Wahrscheinlichkeiten werden definiert.

Zeile 4: Aus den Werten der Zufallsvariable und den zugehörigen Wahrscheinlichkeiten wird ein Dataframe mit 2 Spalten gebildet. Damit man später leichter auf die Spalten zugreifen kann, werden ihnen Namen gegeben. Man kann sogar einen Schritt weitergehen und fordern, dass die Werte und Wahrscheinlichkeiten mit diesem Namen im Dataframe abgelegt werden müssen (siehe weiter unten bei is.randomVariable() ). Man könnte natürlich die Zeilen 1, 2 und 4 zu einer Anweisung zusammenfassen, da man die Objekte values und probs im Folgenden nicht mehr benötigt.

Zeile 5: Die Ausgabe zeigt das Dataframe.

Da die Zufallsvariable coin lediglich zwei Werte besitzt, könnte man sie auch mit einer Binomialverteilung modellieren; aber man kann sich leicht Beispiele vorstellen, die sich nicht auf Verteilungen aus den Basis-Paketen zurückführen lassen.

Überprüfen, ob eine Zufallsvariable vorliegt

In den folgenden Kapiteln werden einige Simulationen und Berechnungen mit Hilfe von Zufallsvariablen durchgeführt, die als Dataframes gespeichert werden. Damit ein Dataframe eine Zufallsvariable beschreiben kann, müssen mehrere Bedingungen erfüllt sein. Es ist daher sinnvoll eine Funktion zu implementieren, die überprüft, ob ein Dataframe diese Bedingungen erfüllt. Die Bedingungen lauten:

Es wird die Konvention vereinbart, dass die Werte und Wahrscheinlichkeiten über die Namen value und prob angesprochen werden. Daher muss das Dataframe zwei Spalten mit diesen Namen besitzen. (Welche die erste und welche die zweite Spalte ist unerheblich, wenn die Spalten nur über die Namen angesprochen werden.)
Die genannten Spalten müssen identische Länge besitzen. Diese Abfrage sollte eigentlich überflüssig sein, da man kein anderes Dataframe erzeugen kann.
Keine der Einzel-Wahrscheinlichkeiten in der Spalte prob darf negativ sein.
Die Summe der Einzel-Wahrscheinlichkeiten muss 1 ergeben. Diese Bedingung kann Schwierigkeiten verursachen, wenn mit den Wahrscheinlichkeiten Berechnungen ausgeführt werden (zum Beispiel wenn Zufallsvariablen miteinander verknüpft werden). Denn durch Rundungsfehler kann es zu kleinen Abweichungen kommen, die man tolerieren kann. Es ist daher besser eine Schranke ε vorzugeben (etwa in der Größenordnung 10^-8) und nur zu überprüfen, ob die Summe der Wahrscheinlichkeiten um weniger als ε von 1 abweicht.
Die Werte in der Spalte value müssen eindeutig sein.

Das folgende Skript zeigt eine mögliche Implementierung für die Funktion is.randomVariable() ; in den cat() -Befehlen kann man geeignete Debug-Ausgaben einfügen oder sie ganz weglassen.

is.randomVariable <- function(rv, eps = 1e-8){
  stopifnot(is.data.frame(rv))
  result <- TRUE
  
  if( !( identical(names(rv), c("prob", "value")) || identical(names(rv), c("value", "prob")) ) ){
    cat("names \n")
    result <- FALSE
  }
  if( !identical(length(rv$value), length(rv$prob)) ){
    cat("length != \n")
    result <- FALSE
  }
  if(any(rv$prob < 0)){
    cat("neg. Prob. \n")
    result <- FALSE
  }
  if( abs( sum(rv$prob) - 1 ) > eps ){
    cat("sum Prob \n")
    result <- FALSE
  }
  if( !identical(rv$value, unique(rv$value)) ){
    cat("Eind. \n")
    result <- FALSE
  }
  return(result)
}

Eingabewerte von is.randomVariable() sind das Dataframe, das eine Zufallsvariable rv (rv = random variable) repräsentieren soll und die Schranke eps.

Man hätte in allen if-Zweigen mit einem return-statement return result die Funktion vorzeitig verlassen können. Der Vorteil der gewählten Implementierung ist, dass alle Verletzungen der Bedingungen an eine Zufallsvariable untersucht und Warnhinweise ausgegeben werden.

Aufgabe: Definieren Sie einige Zufallsvariablen nach dem Vorbild coin oben und testen Sie, ob die Funktion is.randomVariable() den Wert TRUE liefert oder welche der Bedingungen verletzt ist.

Erzeugen von Stichproben mit sample()

Simulationen mit Zufallsvariablen, die durch ein Dataframe gegeben sind

Im Paket base gibt es die Funktion sample(), die verwendet werden kann, um Stichproben zu erzeugen. Sie besitzt 4 Eingabewerte, die sich bis auf replace = FALSE selbst erklären:

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

Der Rückgabewert von sample() ist ein Vektor der Länge size, der aus Elementen von x besteht, die gemäß den Wahrscheinlichkeiten prob "ausgewürfelt" werden.

# Dataframe coin wie oben
sample(x = coin$value, size = 10, replace = TRUE, prob = coin$prob)
# [1] 0 1 1 0 0 1 1 0 0 1

Oder mit einer kleinen Auswertung – man lässt mit Hilfe von table() anzeigen, wie oft die Komponenten von x in der Stichprobe enthalten sind:

s <- sample(x = coin$value, size = 1000, replace = TRUE, prob = coin$prob)
table(s)
# s
# 0   1 
# 377 623

Das folgende Beispiel modelliert einen gezinkten Würfel, bei dem die 6 mit erhöhter und die 1 mit erniedrigter Wahrscheinlichkeit erscheint:

values <- (1:6)
probs <- c(1, 2, 2, 2, 2, 3) / 12

# gezinkter Würfel:
d <- data.frame(value = values, prob = probs)
print(d, digits = 2)
#   value  prob
# 1     1 0.083
# 2     2 0.167
# 3     3 0.167
# 4     4 0.167
# 5     5 0.167
# 6     6 0.250

# Stichproben:
sample(x = d$value, size = 10, replace = TRUE, prob = d$prob)
# [1] 2 1 3 2 5 4 5 6 5 5

s <- sample(x = d$value, size = 1000, replace = TRUE, prob = d$prob)

# absolute Häufigkeiten
table(s)
# s
# 1   2   3   4   5   6 
# 88 158 171 177 138 268

# relative Häufigkeiten
table(s) / 1000
# s
# 1     2     3     4     5     6 
# 0.088 0.158 0.171 0.177 0.138 0.268 

# Mittelwert:
sum(s) / 1000
# [1] 3.923

Die Skripte erwecken womöglich den Eindruck, dass man sample() mit einem Dataframe konfigurieren muss – das ist natürlich nicht der Fall. Man kann die aufrufenden Argumente x und prob auch direkt setzen. Da in vielen Anwendungen x-Werte wie x = (1:6) vorkommen, kann man dies durch x = 6 abkürzen. Das folgende Skript zeigt obige Simulation mit einer Stichprobe der Länge 1000 für einen Laplace-Würfel, der nicht eigens als Dataframe modelliert wird:

s <- sample(x = 6, size = 1000, replace = TRUE, prob = rep(x = 1/6, times = 6))

table(s)
# s
# 1   2   3   4   5   6 
# 162 188 158 181 148 163 

table(s) / 1000
# s
# 1     2     3     4     5     6 
# 0.162 0.188 0.158 0.181 0.148 0.163 

sum(s) / 1000
# [1] 3.454

Man erkennt, dass der Mittelwert jetzt näher bei 3.5 liegt. Abbildung 1 zeigt die Auswertung der Simulationen im Stabdiagramm.

Abbildung 1: Darstellung der mit sample() erzeugten Simulationen (mit jeweils 1000 Würfen) im Histogramm; aufgetragen sind die absoluten Häufigkeiten der Augenzahlen. (Links der gezinkte Würfel, rechts der Laplace-Würfel.)

Wird für x eine natürliche Zahl n vorgegeben, ist der default-Wert für size gleich n:

s <- sample(x = 6, replace = TRUE, prob = rep(x = 1/6, times = 6))
s
# [1] 1 1 3 6 5 5
table(s)
# s
# 1 3 5 6 
# 2 1 2 1

Das Argument replace von sample()

Das Argument replace von sample() ermöglicht

die Simulation von Ziehen ohne Zurücklegen mit replace = FALSE beziehungsweise
Ziehen mit Zurücklegen mit replace = TRUE .

Damit ist Folgendes gemeint:

Ziehen ohne Zurücklegen: Beim Erzeugen einer Stichprobe wird die erste Komponente s₁ aus der Grundmenge x entnommen. Um die zweite Komponente zu ziehen, wird s₁ aus der Grundmenge entfernt und ein Element aus der verkleinerten Grundmenge gezogen. Und so weiter. Aber dann kann die Stichprobe höchstens so viele Komponenten besitzen wie x. Ist size größer als die Länge von x, erhält man eine Fehlermeldung (siehe Beispiel unten).
Ziehen mit Zurücklegen: Jetzt stehen bei jedem Zug alle Elemente von x zur Verfügung; das Argument size kann beliebige Werte annehmen (siehe Beispiele oben).

Dass das Argument replace heißt, sollte man so verstehen: im Falle von replace = TRUE wird nach jedem Zug die Ausgangssituation wiederhergestellt, also die eigentlich verkleinerte Grundmenge x wieder durch das ürsprüngliche x ersetzt.

Beispiel: Ziehen ohne Zurücklegen mit dem gezinkten Würfel

Der Würfel ist mit den Zahlen 1 bis 6 beschriftet. Beim Ziehen ohne Zurücklegen kann eine einmal gezogene Zahl nicht nochmal erscheinen. Das folgende Skript realisiert dies für:

size = 5 : es werden 5 verschiedene Zahlen gezogen (Zeile 4).
size = 6 : die Stichprobe ist eine Permutation von (1, 2, ..., 6) (Zeile 8).
size = 7 : da für den siebten Zug keine Zahl mehr zur Verfügung steht, werden auch die ersten 6 Züge nicht realisiert; es wird eine Fehlermeldung ausgegeben (Zeile 12).

# Dataframe d wie oben

# size = 5:
sample(x = d$value, size = 5, replace = FALSE, prob = d$prob)
# [1] 6 4 2 5 3

# size = 6:
sample(x = d$value, size = 6, replace = FALSE, prob = d$prob)
# [1] 2 3 6 4 1 5

# size = 7:
sample(x = d$value, size = 7, replace = FALSE, prob = d$prob)
# Error in sample.int(length(x), size, replace, prob) : 
#   cannot take a sample larger than the population when 'replace = FALSE'

Beispiel: Ziehung von Lottozahlen

Es soll folgende Simulation des Lottospiels "6 aus 49" durchgeführt werden:

Ein Spieler gibt einen Lotto-Tip ab.
Es werden so lange neue Lottozahlen gezogen bis der Spieler zum ersten Mal 4 Richtige oder mehr als 4 Richtige hat (anstelle der 4 kann man auch jede andere Zahl von 1 bis 6 einsetzen).
Die Folge der erzeugten Lottozahlen soll ausgewertet werden:
- Wie oft wurde gespielt?
- Wie oft wurden "0 Richtige", "1 Richtige", ..., "6 Richtige" erzielt?

Lösungsvorschlag:

Die Funktion numberOfSuccesses() vergleicht einen abgegebenen Lotto-Tip mit den tatsächlich gezogenen Zahlen actual und berechnet die Anzahl der "Richtigen" (Mächtigkeit des Mengendurchschnitts):

# Anzahl der Richtigen 
numberOfSuccesses <- function(tip = (1:6), actual){
  stopifnot( identical(length(tip), length(actual)) )
  return( length( intersect(tip, actual) ) )
}

Die eigentliche Simulation erfolgt in der Funktion simulation(), siehe folgendes Skript. Sie erzeugt so lange neue Lotto-Tips bis die Anzahl der "Richtigen" die Anzahl successes erreicht oder überschreitet. Damit die Simulation reproduzierbar ist, kann man seed für den Zufallsgenerator setzen (Zeile 4). Die eigentlich relevanten Eingabewerte für die Simulation sind:

Die Anzahl der "Richtigen" successes; die Simulation soll abbrechen, wenn diese Zahl erreicht oder überschritten wird (bei successes = 0 muss keine Ziehung stattfinden).
Der abgegebene Lotto-Tip tip.

# so lange wiederholen bis successes 
# oder mehr Treffer erreicht werden; Auswertung
simulation <- function(successes = 4, tip = (1:6), seed = 42){
  set.seed(seed)
  # Anzahl der Treffer im Tip:
  n.succ <- 0
  # Auswertung, Zähler:
  cnt = rep(x = 0, times = 7)

  while(n.succ < successes){
    n.succ <- numberOfSuccesses(tip = tip, actual = sample(x = 49, size = 6))
    cnt[n.succ + 1] <- cnt[n.succ + 1] + 1
    # print(cnt)
  }
  
  # Auwertung:
  cat("Anzahl der Ziehungen: ", sum(cnt), "\n")
  cat("Häufigkeiten der Trefferzahlen:\n")
  print(cnt)
}

Zeile 10 bis 14: Das Kernstück der Simulation ist die while-Schleife. In ihr werden immer wieder neue Lottozahlen gezogen und gezählt, wie viele Richtige der Tip hat. Die Schleife wird wiederholt, solange die Anzahl der Richtigen kleiner ist als successes (Zeile 12).

Zeile 11: Die Ziehung der Lottozahlen wird mit Hilfe von sample() realisiert; dazu werden aus den Zahlen von 1 bis 49 genau 6 Zahlen ausgewählt. Da replace = FALSE (default-Wert) kann eine Zahl nicht mehrfach gezogen werden.

Zeile 8 und 10: Die Variable cnt dient der Auswertung der Simulation. Sie zählt wie oft "0 Richtige", "1 Richtige" und so weiter vorkommen. Dazu wird cnt als Vektor der Länge 7 vorbereitet und alle Komponenten gleich null gesetzt (Zeile 8). In der Schleife wird die entsprechende Komponente um 1 erhöht (Zeile 12); am Ende wird cnt ausgegeben (Zeile 18). Die Gesamtzahl der Ausspielungen erhält man, indem man die Trefferzahlen aus cnt summiert (Zeile 17).

Das folgende Skript zeigt den Aufruf von simulation() und eine typische Auswertung:

simulation(successes = 4)
# Anzahl der Ziehungen:  263 
# Häufigkeiten der Trefferzahlen:
#   [1] 130 102  26   4   1   0   0

Das Argument prob von sample()

Die bisher besprochenen Beispiele werfen einige Fragen auf:

Bei der Modellierung einer Zufallsvariable als Dataframe wurde argumentiert, dass die Summe der Wahrscheinlichkeiten eventuell ein wenig von 1 abweichen kann, insbesondere wenn die Wahrscheinlichkeiten durch vorherige Berechnungen entstanden sind (Fortpflanzung von Rundungsfehlern). Wie reagiert die Funktion sample(), wenn sich die Wahrscheinlichkeiten nicht zu 1 addieren?
Wie werden die neuen Wahrscheinlichkeiten nach Entfernen eines Elementes aus x berechnet, wenn sample() mit replace = FALSE verwendet wird? Es wird sich zeigen, dass dies mit der Antwort auf die erste Frage zusammenhängt.

Im folgenden Beispiel wird die gezinkte Münze von oben mit offensichtlich unsinnigen Wahrscheinlichkeiten definiert und mit ihr eine Simulation durchgeführt. Zum Vergleich wird die identische Simulation mit "korrekt" definierten Wahrscheinlichkeiten durchgeführt:

set.seed(42)
s <- sample(x = c(0, 1), size = 1000, replace = TRUE, prob = c(4, 6))
table(s)
# s
# 0   1 
# 377 623 

set.seed(42)
s <- sample(x = c(0, 1), size = 1000, replace = TRUE, prob = c(0.4, 0.6))
table(s)
# s
# 0   1 
# 377 623

Man erkennt, dass das Argument prob nicht den Anforderungen einer Wahrscheinlichkeit genügen muss. Denn die Zahlenwerte werden als Gewichte interpretiert, die nicht auf 1 normiert sein müssen. Erst wenn die Interpretation als Gewichte versagt, erhält man eine Fehlermeldung:

s <- sample(x = c(0, 1), size = 1000, replace = TRUE, prob = c(0, 0))
# Error in sample.int(length(x), size, replace, prob) : 
#   too few positive probabilities

Aber diese Interpretation erklärt, wie die Wahrscheinlichkeiten neu berechnet werden, wenn beim Ziehen ohne Zurücklegen, also replace = FALSE , ein Element aus x entfernt wird: Es wird einfach das zugehörige Gewicht aus dem Vektor prob entfernt.

Wiederholung von Simulationen mit replicate()

In vielen Anwendungen soll eine Stichprobe nicht nur einmal sondern öfters erzeugt werden; die Auswertung geschieht dann über die Gesamtheit der Stichproben. Naheliegend ist es, die Stichprobe innerhalb einer Schleife zu erzeugen und geeignete Variablen vorzubereiten, die die Ergebnisse der Stichproben abspeichern und die später ausgewertet werden.

Oft kann man die Schleifen durch den geschickten Einsatz der Funktion replicate() vermeiden. Sie wurde in Die Familie der apply-Funktionen in R Teil 3: Weitere mit apply() verwandte Funktionen erläutert und wird hier nicht beschrieben; stattdessen wird ein typisches Beispiel gezeigt, wie man replicate() einsetzen kann.

Beispiel: Zählen der Fixpunkte von Permutationen

Das folgende Listing zeigt die Zahlen von 1 bis 7 sowie 3 Permutationen der Zahlen von 1 bis 7. Eine Permutation eines Vektors enthält alle Komponenten des Vektors, sie können aber in beliebiger Reihenfolge neu angeordnet werden.

1 2 3 4 5 6 7
# Permutationen:
1 2 7 4 6 5 3
4 6 1 5 3 7 2
3 5 2 4 7 1 6

Die erste Permutation enthält 3 Fixpunkte; damit ist gemeint, dass genau 3 Zahlen an der selben Stelle stehen wie im ursprünglichen Vektor (nämlich die Zahlen 1, 2 und 4). Die zweite Permutation enthält keinen Fixpunkt, die dritte Permutation enthält einen Fixpunkt. Und unter dem Blickwinkel der Fixpunkte betrachtet, kann man die erste Zeile des Listings als eine Permutation von (1:7) mit 7 Fixpunkten auffassen.

Es ist ein bekanntes Problem der Kombinatorik zu fragen, wie viele Permutationen von (1:n) es gibt, die genau 0, 1, 2, ..., n Fixpunkte besitzen. Es soll hier nicht dieses Abzählproblem behandelt werden, sondern es soll gezeigt werden, wie man dazu eine einfache Simulation und Auswertung schreibt.

Eine Permutation von (1:n) wird mit Hilfe von sample() erzeugt, wobei man ihre default-Werte verwenden kann:

n <- 7
sample(x = n)
# 7 6 2 5 3 1 4

Um die Anzahl der Fixpunkte einer Permutation zu bestimmen, kann man die Funktion fixedPoints() einsetzen, die sich leicht implementieren lässt:

fixedPoints <- function(p){
  return( sum( p == (1:length(p)) ) )
}

Mit Hilfe von replicate() lasen sich nun beliebig viele Permutationen erzeugen und jeweils die Fixpunkte zählen:

v <- replicate(n = 1000, expr = fixedPoints(sample(x = 7)))

Der Vektor v hat die Länge 1000 und seine Komponenten geben an, wie viele Fixpunkte die entsprechenden Permutationen haben; er lässt sich leicht auswerten:

table(v)
# v
#   0   1   2   3   4   5 
# 359 367 201  58  10   5

Die folgende Abbildung 2 zeigt die Auswertung von 1000 Permutationen von (1:7) beziehungsweise (1:21) .

Abbildung 2: Auswertung von jeweils 1000 Permutationen: gezählt wird, wie viele Fixpunkte sie enthalten. Die relativen Häufigkeiten sind auf der y-Achse aufgetragen.