7 Grundlagen der Simulation in:

Andreas Behr, Ulrich Pötter

Einführung in die Statistik mit R, page 99 - 117

2. Edition 2010, ISBN print: 978-3-8006-3599-3, ISBN online: 978-3-8006-4878-8, https://doi.org/10.15358/9783800648788_99

Series: Vahlens Kurzlehrbücher

Bibliographic information
G S Statistik benötigt Daten. Woher kommen Daten? Eine Möglichkeit, sich Daten zu verschaffen, besteht in der Simulation von Daten. Wir beschreiben in diesem Kapitel die Grundlagen der Simulation, insbesondere die Erzeugung von Folgen von Zufallszahlen und ihre Probleme: Was soll es bedeuten, „Zufall“ auf einem Rechner nachzuahmen? 7.1 Zufallszahlen? 7.2 Gleichverteilung 7.2.1 Monte-Carlo-Integration 7.2.2 Integrale und Erwartungswerte 7.2.3 Endliche Gleichverteilung und Periode 7.2.4 Startwerte 7.3 Zufallszahlen mit vorgegebener Verteilung 7.4 Parametrisierung der Verteilungsklassen 7.5 Stichproben und Tabellen 7.6 Funktionen von Zufallsvariablen 7.7 Übungsaufgaben . Zufallszahlen? Die Begriffe Zufall und Wahrscheinlichkeit sind seit jeher eng mit Glücksspielen verknüpft. Bis heute verweist jede Einführung in die Wahrscheinlichkeitsrechnung auf Münzwürfe, Würfel oder Urnen. Wie simuliert man aber etwaMünzwürfe auf einemComputer? Ein offensichtlicher Unterschied besteht darin, dass Beschreibungen von Münzwürfen viele physikalische Einzelheiten enthalten: die Wurfhöhe undWurfrichtung, die Drehungsgeschwindigkeit der Münze, elastische Effekte der Landung etc. Auch ob das Ergebnis eines Mathematiker und Physiker haben sich vieleGedanken gemacht, warumgerade das physikalisch recht einfache System von Münzwürfen als Prototyp für zufällige Ereignisse geeignet ist. Eine spannende Diskussion, die auch Aspekte von Modellierungsstrategien diskutiert, ist Persi Diaconis: A place for Philosophy?The rise of modelling in statistical science, Quaterly Appl. Math., , , – . 4 10 35 83 _W iS o B eh r P öt te r 2A - B g 4 Grundlagen der Simulation bestimmten Münzwurfs als zufällig erzeugt gelten soll, knüpft an Beschreibungen des Wurf- und Landeverhaltens an. Würfe mit sehr geringer Wurfhöhe, geringer Drehung der Münze oder Landung auf sehr unelastischen Untergründen lassen regelmäßig an der Zufälligkeit der Ergebnisse zweifeln. Aber auf Computern steht nur die Folge von Ergebnissen zur Verfügung. Und man braucht auch nur die Ergebnisse. Kriterien für die Zufälligkeit einer Folge von Ergebnissen können sich also nicht auf das physikalische Vorgehen beim Erzeugen von Ereignissen beziehen. Außerdem verbietet sich der direkte Zugriff auf physikalisch beschreibbare zufällige Ereignisse für Zwecke der Simulation auf Computern, weil dann die Ergebnisse einer Simulation durch andere nicht mehr überprüfbar sind. Wir brauchen also Kriterien, die die verschiedenen Ansprüche für die Simulation an Zufallsfolgen präzisieren. Zunächst ist klar, dass man sich in dieser Perspektive nicht auf Eigenschaften einzelner Zahlen oder Symbole beziehen kann. Vielmehr muss man sich auf Folgen von solchen Zahlen beziehen. Sie sollen zudem auf dem Rechner dargestellt werden. Folglich ist die Folge immer endlich, weil man nur eine endliche Anzahl von Zahlen auf dem Rechner darstellen kann. Ein Kriterium, dass man an „zufällige“, auf Rechnern realisierte Zahlenfolgen stellen sollte, wurde schon genannt: Sie sollten sich reproduzieren lassen. Man sollte in der Lage sein, seine Ergebnisse auch nach einer Woche oder einem Jahr reproduzieren zu können und andere (auf anderen Rechnern) sollten nach der Beschreibung des Programms in der Lage sein, exakt gleiche Ergebnisse zu erhalten. Das ist der entscheidende Unterschied zu Ergebnissen von Glücksspielen, die sich eben nicht wiederholen lassen. Folglich braucht man Algorithmen, rechnerische Verfahren, die die Zahlenfolgen erzeugen. Damit ist auch ein weiteres Kriterium zufälliger Folgen ausgeschlossen, das in der philosophischen Diskussion oft angeführt wird: Die Folgen, die wir benutzen wollen, sind sicherlich nicht überraschend oder unvorhersagbar. Kennt man den Algorithmus sowie die Startwerte, dann kennt man alle weiteren Werte der erzeugten Zahlenfolge. Ein zweites Erfordernis scheint ebenfalls klar zu sein: Die erzeugten Zahlenfolge sollte eine angebbare Verteilung haben. Hat man einen Algorithmus, der die Folge x , x , . . . erzeugt, dann sollte lim n n n i= 1(a xi b) = F(b) F(a) (7.1) Kriterien der „korrekten“ Durchführung von Münzwürfen, für „korrektes“ Vorgehen beim Würfeln, beim Roulette, beim Mischen von Spielkarten oder bei Lotto-Ziehungen beziehen sich immer auf das Verfahren, mit dem gewürfelt, gemischt oder gezogen wird. Auch wenn Einführungen in die Statistik gern anderes nahelegen, werden statistische Tests als Kriterien praktisch nicht genutzt. Denn statistische Tests können sich nur auf die Ergebnisse der Münzwürfe etc. beziehen. Daher erlauben sie nur indirekt Rückschlüsse auf das Verfahren. Es gibt natürlich auch Verfahren der Erzeugung von zufälligen Folgen, die an dem zufälligen Charakter von physikalischen Ereignissen ansetzen. Das Paket random stellt etwa einen automatisierten Zugriff auf Zahlenfolgen bereit, die durch atmosphärisches Rauschen erzeugt werden. Solche zufälligen Folgen haben ihre Anwendung in der Kryptographie und ähnlichen Bereichen, sind aber für unsere Überlegungen weitgehend irrelevant. . Gleichverteilung für eine vorgegebene Verteilungsfunktion F(.) gelten. Wir betrachten zunächst Gleichverteilungen auf dem Intervall [ , ]. Dann ist der obige Grenzwert gleich b a, wenn a, b [ , ] sind. Ein drittes Erfordernis ist sicherlich, dass der rechnerische Aufwand zur Erzeugung einer zufälligen Folge von Zahlen sehr gering ist, denn man benötigt oft sehr viele Zufallszahlen, etwa wenn man einige Millionen Datensätze mit jeweils tausenden Beobachtungen erzeugen möchte, um statistische Verfahren an verschiedenen Datensätzen auszuprobieren. . Gleichverteilung Nun ist es erstaunlich, dass relativ einfache mathematische Funktionen diese drei Kriterien erfüllen. Setzt man z.B. ui := i a i a (7.2) dann ist die Folge ui gleichverteilt auf dem Intervall [ , ], wenn die Konstante a irrational ist. Ein entsprechendes Verfahren lässt sich näherungsweise in R umsetzen. Zwar kann man keine irrationalen Zahlen in R (oder einem anderen Computerprogramm) exakt numerisch darstellen, aber man kann Näherungen ausprobieren. Nimmt man etwa a = π, dann würde eine Näherung in R so aussehen: Das Symbol 1(A) steht für die Indikatorvariable des Ereignisses A, d.h. sie nimmt den Wert an, wenn A wahr ist, sonst. Für endliche Folgen kann diese Forderung nur approximativ gelten. Aber wir werden sehen, dass einige praktisch verwandte Verfahren das Kriterium erfüllen, wenn unbegrenzter Speicherplatz für die Zahlendarstellung unterstellt wird. Generell lassen sich auch alle (berechenbaren) statistischen Tests als Kriterien für die „Zufälligkeit“ von Zahlenfolgen heranziehen. Das führt zu einem möglichen Kriterium für die mathematische Definition von „zufälligen“ Folgen. Überblicke über Versuche, „Zufälligkeit“ zumindest für Folgen von Zahlen zu definieren, geben Sergio Volchan:What is a random sequence? Am. Math. Monthly, , January , – und Andrei A. Muchnik et al.: Mathematical metaphysics of randomness,Theoretical Computer Science, , , – . Sogar die Dezimalentwicklung (oder die Entwicklung in anderen Basen) etlicher mathematischer Konstanten scheint diese drei Kriterien zu erfüllen. Empirisch ist die Anzahl der Ziffern - in der Dezimaldarstellung vieler Konstanten gleichverteilt. Seit Bailey, Borwein und Plouffe eine Methode entwickelt haben, die Darstellung einer beliebigen Stelle der Entwicklung von π und anderen Konstanten anzugeben, ohne die vorhergehenden zu berechnen, scheint der Weg offen, statistische Eigenschaften dieser Konstanten zu erforschen. Einen kurzen Überblick geben David Bailey und Jonathan Borwein: Experimental mathematics and computational statistics, WIREs Computational Statistics, vol. , , – . u ist die größte ganze Zahl kleiner als u. Die ui sind also die Zahlen hinter dem Dezimalkomma von i a. Ist etwa a = π = . ..., die Kreiszahl, dann ist u = . .... Dass diese Folge gleichverteilt ist, wurde von HermannWeyl bewiesen. Einen Überblick über ähnliche Verfahren und verwandte Ergebnisse werden im Buch von Walter Freiberger und Ulf Grenander: A Short Course in Computational Probability and Statistics, Springer, , vorgestellt. Es gibt eine enge Verbindung mit der Frage der Gleichverteilung von Ziffern in mathematischen Konstanten (David Bailey und Richard Crandall: On the random character of fundamental constant expansions, Experimental Mathematics, , , – ). Grundlagen der Simulation > rng < function(n){neu < (1:n∗pi) + neu Ćoor(neu)} > rng(5) [1] 0.1415927 0.2831853 0.4247780 0.5663706 0.7079633 Folgen von Zahlen, die einige Kriterien von Zufallsfolgen erfüllen und auf dem Rechner dargestellt werden können, nennt man Pseudozufallszahlen. . . Monte-Carlo-Integration Die Eigenschaft der Gleichverteilung (Kriterium mit F(u) = u) reicht nun schon für eine wichtige Aufgabe der Simulation: Die numerische Berechnung von Integralen. Das Problem, Werte von Integralen zu berechnen, entsteht etwa, wenn man etwas über das durchschnittliche Verhalten von zufälligen Systemen sagen möchte. Im einfachsten Fall wird das System nur durch eine interessierende Größe charakterisiert, etwa die Wartezeiten bis zu einem bestimmten Ereignis, wie das Eintreffen von Bussen oder Bahnen, dieWartezeit an den Kassen von Supermärkten oder in Computernetzwerken oder auch die Dauer von Ehen. Die klassische mathematische Methode, solche Größen zu beschreiben, besteht darin, sie als Zufallsvariable aufzufassen, als eine Funktion, die jedem ω Ω eine entsprechende Dauer zuweist. Formal würde man für solche Wartezeiten also schreiben: X : Ω [ , ) Nun ist die genaue Form von Ω weitgehend irrelevant, soweit die Menge nur hinreichend viele Elemente hat. Wichtig ist nur, dass sie zusätzlich eine Wahrscheinlichkeitsstruktur besitzt, eine weitere Abbildung, die für eine hinreichend große Klasse von Teilmengen von Ω diesen Teilmengen Wahrscheinlichkeiten zuordnet. Für einfache Systeme kann man daher Ω = [ , ] wählen und als Wahrscheinlichkeitsmaß die Gleichverteilung auf [ , ]. Die durchschnittliche Wartezeit, der Erwartungswert der Zufallsvariablen X, ist dann der Durchschnitt aller Wartezeiten, gewichtet mit ihren „Wahrscheinlichkeiten“. Da man die Wartezeiten nicht aufzählen kann, also einfache Durchschnitte nicht berechnen kann, ersetzt man die Durchschnittsbildung (endliche gewichtete Summen) durch passende Integrale. Formal schreibt man E(X) := Ω X(ω) Pr(dω) Jede reelle Zufallsvariable lässt sich so darstellen. . Gleichverteilung Hat man aber Ω = [ , ] und als Wahrscheinlichkeitsmaß die Gleichverteilung gewählt, dann wird daraus E(X) = X(u) du Hätte man eine zufällige Folge U ,U , . . . mit Werten in [ , ], die den Gesetzen der Wahrscheinlichkeitstheorie folgt, dann würde das starke Gesetz der großen Zahlen implizieren, dass lim n n n i= X(Ui) = E(X) mit Wahrscheinlichkeit gelten würde. Für ein hinreichend großes n könnte man also /n X(ui) als Näherung für E(X) betrachten. Die übrigen Gesetze der Wahrscheinlichkeitstheorie wie der zentrale Grenzwertsatz und der Satz vom iterierten Logarithmus erlauben dann die Abschätzung des Näherungsfehlers. Nun zeigt sich, dass die Wahl der Folge ui = i a i a zumindest ebenso gut zur Näherung des Integrals E(X) geeignet ist wie „echte“ Zufallszahlen. Denn der Grenzwert der Durchschnitte derX(ui) in der obigen Formel gilt für jede Folge ui der Form (7.2), solange nur a irrational ist. Insbesondere gilt die Konvergenz nicht nur mit Wahrscheinlichkeit . Es ist also gar kein Verlust, dass unsere Folge sicherlich keine „echte“ Zufallsfolge ist. Methoden der numerischen Berechnung von Integralen mit Hilfe von Pseudozufallszahlen werden Monte-Carlo-Methoden genannt. . . Integrale und Erwartungswerte Wir probieren das Verfahren einmal aus und überlegen uns ein einfaches Modell für Wartezeiten. Wartezeiten sind sicherlich positiv, also könnte die Verteilungsfunktion vonX etwa F(x) = exp( λx) wie in Abbildung . sein. Diese Verteilungsfamilie nennt man Exponentialverteilung. Wir müssen zunächst sehen, wie die AbbildungX : Ω [ , ) aussehen soll, damitX die Verteilungsfunktion F(.) hat. Wenn U gleichverteilt auf [ , ] ist und fürX(.) die Quantilsfunktion F (u) wie in Abschnitt . . gewählt Das folgt leicht aus dem Beweis der Gleichverteilung im schon erwähnten Buch von Freiberger und Grenander. Mit etwas mehr technischem Aufwand zeigen sie auch, dass der Fehler bei der Näherung schneller als /n gegen geht, also sogar schneller als für jede „echte“ Folge unabhängiger Zufallszahlen, wenn auch nur für Irrationalzahlen amit bestimmten Eigenschaften. Die guten Eigenschaften solcher Folgen für die Berechnung von Integralen haben dazu geführt, etliche „deterministische“ Folgen mit ebenso guten oder besseren Eigenschaften als Folgen von Zufallszahlen zu entwickeln, die auch allgemeine Integrationsprobleme behandeln können. Einen kurzen Überblick liefert z.B. Rüdiger Seydel: Tools for Computational Finance, Springer , Kapitel . Grundlagen der Simulation wird, erhält man das gewünschte Ergebnis. Für streng monotone F(.) wie die Exponentialverteilung sieht man das auch gleich, denn Pr(X(U) x) = Pr(F (U) x) = Pr(U F(x)) = F(x) Abbildung 7.1: Exponentialverteilungen. Verteilungsfunktion für λ = und λ = (gepunktet). Die letzte Gleichung folgt aus der Definition der Gleichverteilung, für die Pr(U z) = z für z [ , ] ist. Will man also den Erwartungswert vonX ausrechnen, dann muss man das Integral E(X) = F (u) du = log( u)/λ du berechnen. Dieses Integral lässt sich natürlich relativ leicht analytisch berechnen und das Ergebnis ist /λ. Wir benutzen nundie Folgeui = i π i π , um das Ergebnis numerisch zu berechnen. Wir benutzen die vorher definierte Funktion rng() und berechnen denDurchschnitt der FunktionswerteM( X(ui) i = , . . . , n ), etwa für n = , , , . Das Ergebnis ist > exponential < function(u,lambda=1){ log(1 u)/lambda} > f10 < mean(exponential(rng(10),lambda=1));f10 [1] 1.072721 > f100 < mean(exponential(rng(100),lambda=1));f100 [1] 0.9942059 > f1000 < mean(exponential(rng(1000),lambda=1));f1000 [1] 1.045532 > f1000000 < mean(exponential(rng(1000000),lambda=1));f1000000 [1] 1.000171 In diesem Beispiel erhält man Ergebnisse mit einer Genauigkeit von - Stellen mit recht geringem Rechenaufwand. Für eine Verbesserung der Näherung ist aber ein erheblich höherer Rechenaufwand nötig. Der Vorteil der Monte-Carlo-Methode der Integration ist offensichtlich: Sie kann ohne weiteres für beliebig komplizierte Funktionen benutzt werden, denn sie erfordert keine Anpassung an die Form der zu integrierenden Funktion. Die Nachteile sind auch klar: Die Konvergenz kann sehr langsam sein, und selbst für relativ geringe Genauigkeit ist u.U. ein erheblicher Rechenaufwand notwendig. Zudem kann die Fehlerabschätzung in den meisten Fällen nur mit statistischen Methoden erreicht werden, ist also erstens selbst unsicher und beruht zweitens nur auf einer Analogie zu „echten“ Zufallszahlen, die für die verwandten Pseudozufallszahlen nicht zutreffen müssen. Beide Schwierigkeiten lassen sich zwar oft umgehen, erfordern dann aber wieder eine genauere Analyse entweder der zu integrierenden Funktion oder der verwandten Pseudozufallszahlen. . Gleichverteilung . . Endliche Gleichverteilung und Periode Abbildung 7.2: Aufeinanderfolgende Pseudozufallszahlen nach Formel (7.2). Ein Problem der Folge ui lässt sich leicht illustrieren: Im Gegensatz zu Folgen vonWerten unabhängiger Zufallsvariabler sind die aufeinanderfolgendenWerte der Folge ui sicher nicht unabhängig: Malt man ein Bild der aufeinanderfolgenden Werte, etwa mit > a < rng(50) > gerade < 1:25∗2 > ungerade < gerade 1 > plot(a[ungerade],a[gerade],type="p") dann erhält man Abbildung . . Es ist klar, dass die nächste erzeugte Zahl ui+ gerade ui + . . . . ist, solange ui < . . . . gilt. Ist ui . . . . hat man ui+ = ui + . . . . . Solche Folgen sind zwar gut geeignet, einfache Integrale zu berechnen, aber sie sind offenbar völlig ungeeignet, Zeitreihen zu simulieren oder Integrale von Funktionen mehrerer Argumente zu berechnen. Denn zwei aufeinanderfolgende Zahlen der Folge ui sind stark korreliert. Für die Berechnung von Integralen in höheren Dimensionen bedeutet das, dass nur ein kleiner Teil des Definitionsbereichs der zu integrierenden Funktion durchlaufen wird. Untersucht man etwa Funktionen, die von Werten im Quadrat wie in Abbildung . abhängen, dann ist klar, dass allein die Funktionswerte entlang der beiden Geradenstücke berechnet werden. Man benötigt also für solche Probleme eine zusätzliche Eigenschaft von Pseudozufallszahlen, die eine gleiche Verteilung von Punkten auch in , , . . . Dimensionen sicherstellen. Eine wichtige Kenngröße für Pseudozufallszahlen ist daher die größte Dimension, bis zu der eine Gleichverteilung erreicht werden kann. Bisher haben wir so getan, als könne man eine irrationale Zahl, etwa π, exakt auf einem Rechner darstellen. Das geht aber nicht, weshalb Implementationen von Folgen ui wie (7.2) durch die Funktion rng() nicht alle Eigenschaften haben können, die der Folge ui zukommen. Die Gleichverteilung der Folge hängt daran, dass nach jeder aufsteigenden Teilfolge die nächste um ein wenig versetzt ist. Die Größe dieser Versetzung hängt von den versetzten Dezimalstellen von π ab. Da π irrational ist, wiederholen sich diese Zahlen nie, man beginnt immer Formal betrachtet man k-fache Versionen für den Grenzwert in (7.1). Dazu ersetzt man das Intervall (a, b) durch k Intervalle (ai , bi) und betrachtet die Häufigkeit der Ereignisse a < xi b , . . . , ak < xi+k bk . Der Grenzwert soll dann dem Produkt über alle Terme F(bi) F(ai) entsprechen. Anschaulich sollte der Anteil der Folgenglieder, die in ein Rechteck mit den Koordinaten (ai , bi) fallen, dessen Volumen entsprechen. Damit möchte man entsprechende einfache Unabhängigkeitsforderungen zwischen Elementen der Folge sicherstellen (vgl. etwa Donald Knuth:The Art of Computer Programming, vol. : Seminumerical Algorithms, Addison-Wesley , S. ff). Grundlagen der Simulation an verschiedenen Stellen. Kann man aber nur endlich viele Zahlen zwischen und darstellen, dann muss sich die Stelle irgendwann wiederholen. Damit aber beginnt die ganze Folge von neuem, es wird ja jedesmal ein fixer Betrag (in der Rechnergenauigkeit) addiert. Für Generatoren, bei denen Folgenwerte nur von dem Vorwert (und Konstanten) abhängt, ist die maximale Länge einer Folge, bevor sie sich wiederholt (ihre Periode) gerade die maximale Anzahl der darstellbaren Zahlen. Ist aber die Periode relativ klein, dann können auch nur relativ kurze Folgen von Zufallszahlen verwandt werden. Ein letztes hier betrachtetes Kriterium für möglichst allgemein verwendbare Pseudozufallszahlen ist daher eine möglichst lange Periode. runif Abbildung 7.3: Aufeinanderfolgende Pseudozufallszahlen mit runif(). Bis in die er Jahre wurden fast ausschließlich Generatoren verwandt, die ähnlich wie die Folge (7.2) nur von den Vorwerten abhängen. Da man auch für möglichst kurze Rechenzeiten sorgen wollte, benutzte man die Wortlänge der Prozessoren, also i.d.R. Bit. Damit war deren Periode maximal . . . , oft aber auch viel kleiner. Ein weiteres Problem dieser Generatoren ist, dass je zwei aufeinanderfolgende Werte der Folge verschieden sein müssen. Wären sie es nicht, bliebe die Folge immer beim gleichenWert. Bei einer zufälligen Folge mit diskreten Werten ist das unmöglich, denn in ihnen muss es beliebig lange Teilfolgen gleicher Werte geben. Rs Befehl runif() benutzt in der Voreinstellung einen Mersenne Twister genannten Generator, dessen Periodenlänge ist. Weil die Periode um ein vielfaches länger ist, als es darstellbare Zahlen gibt, entstehen auch Folgen gleicher Zahlen. Die Gleichverteilung ist bis zu Dimensionen garantiert. Bei Dimensionen ergeben sich Abhängigkeiten zwischen diesen Werten. Für Werte in zwei Dimensionen erhält man Abbildung . durch: > a < runif(100) > gerade < 1:50∗2 > ungerade < gerade 1 Auch die besseren in R implementierten Generatoren produzieren z.B. eine zu kleine Zahl von Dubletten im Vergleich zu ihrer wahrscheinlichkeitstheoretisch berechneten Anzahl, wenn gleichverteilte unabhängige Zahlen bis zur Rechnergenauigkeit gerundet werden. Eine Beispielrechnung findet sich in der Hilfe des Befehls set.seed(). Das ist eine riesige Zahl. Die Zahl der Partikel im Universum wird auf nur geschätzt. Allerdings kommen solch enorme Zahlen manchmal in der Statistik vor. So ist die Zahl der verschiedenen Stichproben vom Umfang aus Millionen Bundesbürgern etwa . . Gleichverteilung > plot(a[ungerade],a[gerade],type="p") Der Befehl runif() nimmt als erstes Argument die Anzahl der zu erzeugenden Pseudozufallszahlen. Zwei weitere Argumente können benutzt werden, um ein anderes Intervall [a, b] als [ , ] zu wählen. Andere Generatoren Der Mersenne Twister dürfte für die meisten Anwendungen vollständig ausreichend sein. Für einige spezielle Anwendungen kann es aber notwendig sein, andere Generatoren zu wählen. Im Basispaket base sind fünf weitere Generatoren verfügbar, die durch die Funktion RNGkind() gewählt werden können. Diese Generatoren werden im wesentlichen bereitgestellt, um alte Ergebnisse reproduzieren zu können. Varianten der Generatoren aus alten R-Versionen können ebenfalls gewählt werden. Diese Versionen haben aber alle schlechtere Eigenschaften als der Mersenne Twister. Auf der anderen Seite stellen mehrere Pakete spezielle Generatoren mit interessanten und hilfreichen Eigenschaften zur Verfügung. Das Paket SuppDist stellt den Generator rMWC1019 mit einer Periode von und einer Gleichverteilung auf Dimensionen bereit. Der Generator kann etwa dann notwendig sein, wenn man sehr lange multivariate Zeitreihen simulieren möchte, bei denen die Dimension des Mersenne Twisters ein Problem sein kann. In Extremfällen kann man den Generator randaes des gleichnamigen Pakets verwenden. Dieser Generator erzielt eine Gleichverteilung auf etwa Dimensionen. Zudem werden nicht Bit Wörter erzeugt, sondern Bits, die volle Länge der Mantisse in -Byte Zahlen, dem Standarddatentyp in R. Nur benötigt dieser Generator fast doppelt so viel Rechenzeit wie der Mersenne Twister. Man wird diesen Generator wohl nur verwenden, wenn man besonders erstaunliche Simulationsergebnisse aufklären möchte. Eine wichtige Anwendung, die die Wahl eines alternativen Generators erzwingt, ist die Benutzung mehrerer Prozessoren. Die meisten heutigen Rechner haben mehrere Prozessoren und gerade bei Simulationen ist es naheliegend, die Berechnungen parallel auf die vorhandenen Prozessoren zu verteilen. Die bisher vorgestellten Generatoren produzieren aber einen einzigen Strom von Pseudozufallszahlen. Würde man naiverweise die gleichen Generatorenbefehle in den Teilprozessen benutzen, dann erhielten alle Teilprozesse die gleiche Folge von Pseudozufallszahlen. Man kann also nur zunächst auf einem Prozessor alle notwendigen Pseudozufallszahlen erzeugen und sie dann auf die Teilprozesse verteilen. Das ist nicht nur zeitaufwendig und fehleranfällig, die Aufteilung mag auch zu Problemen mit der Gleichverteilung führen, die zudem von den Einzelheiten der Implementation der Generatoren abhängen. Das Paket rstream stellt den Generator MRG32k3a bereit. Der Generator hat zwar „nur“ eine Periode von und ist gleichverteilt auf mindestens Dimensionen. Aber die Zufallsfolge lässt sich sehr schnell auch in Blöcken berechnen Grundlagen der Simulation und man kann relativ einfach „unabhängige“ Ströme von Pseudozufallszahlen erzeugen. Das Paket rstream enthält auch Methoden für den einfachen Zugriff auf die verschiedenen Ströme sowie für die Aufteilung der Ströme auf mehrere Prozessoren. . . Startwerte Es fehlt noch ein Mechanismus, die Reproduzierbarkeit der Prseudozufallszahlen zu garantieren. Der Befehl runif() liefert in einer Folge von Anwendungen jeweils die nächsten Zahlen eines Stroms von Pseudozufallszahlen, die beim Start von R entweder durch die Uhr des Rechners oder durch den gespeicherten Wert des letzten Zustands des Generators initialisiert werden. Folglich ergeben mehrere Aufrufe des Generators verschiedenWerte und damit werden sich auch die Ergebnisse von Simulationen unterscheiden. Diese Art der Variabilität von Simulationsergebnissen erschwert aber die Fehlersuche in Programmen, die Reproduktion eigener Ergebnisse nach einiger Zeit und die Nachvollziehbarkeit für Andere. Simulationsergebnisse wären für wissenschaftliche Zwecke praktisch wertlos. Man kann die Startwerte der Generatoren in R vorgeben undwürde so die Reproduzierbarkeit der Ergebnisse garantieren können. Nur erfordern verschiedene Generatoren verschiedene Anzahlen von Startwerten (der Mersenne Twister hat Startparameter, der Generator MRG32k3a ), zudem liefern nicht alle Kombinationen der Parameter gute Pseudozufallszahlen. Deshalb stellt R einen einfachenMechanismus zur Verfügung, die Startwerte von Generatoren einheitlich zu wählen. Der Befehl set.seed() hat als erstes Argument eine ganze Zahl, die je nach verwandten Generator die Startwerte festlegt. Ein fix gewählter Wert führt bei gegebenem Generator immer zu gleichen Folgen von Pseudozufallszahlen: > set.seed(3543) > runif(3) [1] 0.8545323 0.2785417 0.2434567 > set.seed(3543) > runif(3) [1] 0.8545323 0.2785417 0.2434567 Wenn man also Simulationsresultate kommunizieren will, dann sollte man immer zunächst den set.seed() Befehl verwenden und ihn auch dokumentieren. Versuche mit verschiedenen Startwerten in set.seed() sind natürlich auch Eine Beschreibung einschließlich konkreter Beispiele findet sich in Pierre L’Ecuyer und Josef Leydold: rstream: Streams of random numbers for stochastic simulation, R News, , , – . Man kann den letzten Zustand des Generators speichern. Er ist in der Variablen .Random.seed abgelegt. Daher würde er bei der nächsten Sitzung mit allen anderen Variablen wieder geladen, wenn man beim Beenden der letzten Sitzung den „Workspace“ gespeichert hat. Das ist allerdings nur hilfreich, wenn man Simulationen zunächst unterbrechen und später (ohne weitere Rechnungen) fortsetzen möchte. . Zufallszahlen mit vorgegebener Verteilung hilfreich, um die Abhängigkeit nicht nur von den Startwerten sondern auch von den Generatoren zu minimieren. . Zufallszahlen mit vorgegebener Verteilung Wir haben zu Beginn dieses Abschnitts schon eine allgemeine Möglichkeit angedeutet, aus gleichverteilten Pseudozufallszahlen Zufallszahlen beliebiger Verteilungen zu generieren. Ist nämlich die Verteilungsfunktion der zu erzeugenden Zufallszahlen F(.), dann hat F (U) diese Verteilung, wenn U gleichverteilt ist. Man kann also beliebig verteilte Pseudozufallszahlen erzeugen, wenn ein Generator für gleichverteilte Pseudozufallszahlen zur Verfügung steht. Nur ist die Quantilsfunktion F (.) oft nicht einfach anzugeben. Zudem gibt es oft effizientere Methoden, Pseudozufallszahlen vorgegebener Verteilungen zu produzieren. Daher stellt R für viele bekanntere Verteilungen spezielle Befehle zur Verfügung. Zudem gibt es eine einheitliche Benennung der entsprechenden Befehle. Zu jeder Verteilungsfamilie gibt es immer Funktionen, die die Dichte, die Verteilungsfunktion, die Quantilsfunktion und entsprechend verteilte Pseudozufallsfolgen erzeugen. Die entsprechenden Befehle beginnen mit d für die Dichten f (.), p für die Verteilungsfunktion F(.), q für die Quantilsfunktion F (.) und r für den entsprechenden Pseudozufallszahlengenerator. So erhält man das berühmte . % Quantil der Standardnormalverteilung, das bei statistischen Tests auf dem % Niveau immer zitiert wird, durch > qnorm(0.975) [1] 1.959964 Wir können das mit der entsprechenden Verteilungsfunktion überprüfen, denn es sollte F(F ( . )) = . gelten: > pnorm(qnorm(0.975)) [1] 0.975 Und wir erhalten normalverteilte Pseudozufallszahlen durch rnorm() > set.seed(22) > rnorm(3) [1] 0.5121391 2.4851837 1.0078262 Das Basispaket stats enthält z.Z. die Funktionen für folgende Verteilungen: Die Beta, Binomial, Cauchy, χ , Exponential, F-, Γ- und Normalverteilungen sowie die geometrische, hypergeometrische, lognormale, logistische, multinomiale, negativ binomiale, Poisson, Student t, Weibull und Wilcoxon-Verteilung. Die entsprechenden Namen der Zufallszahlengeneratoren sind also rbeta(), rbinom() etc. Eine Reihe von Generatoren für weitere Verteilungsfamilien (zusammen mit Dichte-, Verteilungs- und Quantilsfunktionen) sind im Paket SuppDist Grundlagen der Simulation enthalten. Viele andere Pakete enthalten entsprechende Funktionen für spezielle Verteilungen, die in bestimmten Anwendungsfeldern von Bedeutung sind. Pseudozufallszahlen vorgegebenerVerteilung erleichtern die Berechnung von Erwartungswerten und anderen Integralen. Denn an Stelle des Durchschnitts über die Funktion exponential(), der Quantilsfunktion der Exponentialverteilung, zusammen mit einer Folge gleichverteilter Variabler kann nun einfach der Durchschnitt über eine Folge von exponentialverteilten Zufallszahlen berechnet werden. Betrachten wir noch einmal den Erwartungswert einer Exponentialverteilung aus dem letzten Abschnitt. Sie hat die Verteilungsfunktion F(x) = exp( λx). Pseudozufallszahlen können wir nun direkt mit dem Befehl rexp() simulieren. Für den Erwartungswert gilt: E(X) = x dF(x) n n i= xi wobei die xi eine Folge exponentialverteilter Pseudozufallsvariabler mit Ratenparameter λ sind. Die Erwartungswerte lassen sich also einfach durch den Durchschnitt über entsprechend verteilte Pseudozufallsfolgen approximieren. Für λ = kann man also einfach schreiben: > set.seed(88) > mean(rexp(10)) [1] 1.560786 > mean(rexp(100)) [1] 0.8867868 > mean(rexp(1000)) [1] 1.035351 > mean(rexp(1000000)) [1] 0.999259 Für kleine n ist diese Methode offenbar deutlich schlechter als die Methode mit der deterministischen Folge (7.2). Auf der anderen Seite kann man direkt eine Simulationsmethode für den Erwartungswert von beliebigen Funktionen der Zufallsvariablen programmieren, ohne zunächst die Funktionen in äquivalente Ausdrücke mit gleichverteilten Argumenten übersetzen zu müssen. . Parametrisierung der Verteilungsklassen Die Arbeit mit den verschiedenen vordefinierten Zufallszahlengeneratoren erfordert ein minimales Verständnis der jeweils fest vorgegebenen Parametrisierung der Verteilungen. Die Parametrisierungen sind in den entsprechenden Hilfeseiten dokumentiert, sie müssen aber nicht mit den Versionen übereinstimmen, die in anderen Quellen und Lehrbüchern angegeben werden. Die Umrechnung von einer Parameterdarstellung in eine andere ist in den meis- . Stichproben und Tabellen ten Fällen sehr einfach, erfordert aber eine gewisse Aufmerksamkeit für die dokumentierte Variante, die in dem benutzten Paket implementiert ist. Z.B. benutzt der rexp() Befehl die oben bereits angegebene Parametrisierung Pr(X > x) = exp( λx), verschiedene Werte von λ können also durch rexp(n,rate=λ) angegeben werden. In vielen Anwendungen benutzt man aber die Form F(x) = exp( x/µ), denn in dieser Form ist der Parameter µ der Erwartungswert der Exponentialverteilung. Die Umrechnung ist aber denkbar einfach: λ = /µ. . Stichproben und Tabellen Die Simulation von Tabellen kann zwar oft durch die Funktionen für die Simulation multinomialer oder Poisson-Verteilungen erreicht werden. Diese Möglichkeit ist aber wenig benutzerfreundlich, weil der gesamte Vektor der Wahrscheinlichkeiten angegeben werden muss. Die Auswahl von Teilmengen aus endlichen Mengen wird durch diese Funktionen gar nicht repräsentiert. Für solche Fragen stellt R den Befehl sample() zur Verfügung. Der Befehl > set.seed(7) > sample(1:7,3) [1] 7 3 1 zieht aus den Zahlen , . . . , drei Zahlen ohne Zurücklegen. Der Befehl simuliert also das Ziehen von Kugeln aus Urnen, wenn Kugeln nicht zurückgelegt werden. Das entspricht etwa den Lottoziehungen in Deutschland oder der Ziehung von Befragten aus Registern von Bewohnern eines Landes. Das erste Argument des sample() Befehls gibt den Bereich an, aus dem Stichproben gezogen werden sollen. Ist nur eine ganze Zahl n angegeben, dann wird 1:n als Bereich benutzt. Es kann aber auch ein beliebiger Vektor vorgegeben werden. Dann werden Werte aus diesem Vektor gezogen. Das Argument replace mit den Werten TRUE bzw. FALSE gibt an, ob mit oder ohne Zurücklegen der Kugeln gezogen werden soll. Man kann auch Ziehungswahrscheinlichkeiten für einzelne Elemente der Grundmenge durch die Option prob festlegen. Zwar kann man mit dem sample() Befehl auch diverse Tabellen höherer Dimension simulieren. Wenn man aber Tabellen mit vorgegebenen marginalen Verteilungen simulieren möchte, kann der Befehl sample() nicht benutzt Wenn der Befehl sample() innerhalb einer Simulation verwandt wird, muss man darauf achten, dass der Befehl zweideutig werden kann. Möchte man aus einer bereits zufällig gewählten Menge wiederum ein Element auswählen, dann kann es passieren, dass die erste Menge nur noch ein Element n enthält. Ist dieses Element durch eine ganze Zahl repräsentiert, dann wird nicht das einzig vorhandene Element ausgewählt. Stattdessen wird eine Stichprobe aus : n Zahlen gezogen. Das kann zu schwer erkennbaren und „zufälligen“ Fehlermeldungen im weiteren Programm führen. Das Paket sampling enthält etliche weitere Routinen, die häufig verwandte Stichprobendesigns simulieren. Grundlagen der Simulation werden. In diesem Fall kannman aber den Befehl r2dsample() benutzen. Solche Probleme tauchen etwa auf, wenn man Wählerwanderungen zwischen Wahlen simulieren möchte. Dann sind i.d.R. die jeweiligen Wahlergebnisse bekannt undmanmöchte etwas über Veränderungen derWahlentscheidungen zwischen den Wahlen erfahren. Während die Randverteilungen (die Ergebnisse der beiden Wahlen) genau bekannt sind, gilt das nicht für die Wählerwanderungen. Also würde man gern die Ergebnisse der beiden Wahlen konstant halten, aber die möglichen Übergänge simulieren. Das erste Argument des r2dtable ist die Anzahl der zu erzeugenden Tabellen. Die nächsten beiden sind jeweils Vektoren mit den Randverteilungen. Das Ergebnis ist eine Liste mit Matrizen als Elementen. > r2dtable(1,c(40,60),c(50,50)) [[1]] [,1] [,2] [1,] 17 23 [2,] 33 27 . Funktionen von Zufallsvariablen In vielen Fällen muss man statt mit einer Zufallsvariablen mit Funktionen einer oder mehrerer Zufallsvariablen arbeiten. Wenn man etwa wieder an Wartezeiten denkt, dann interessiert oft die Gesamtwartezeit, die sich aus der Summe der Wartezeiten an einzelnen Stationen ergibt. Man möchte also über die Zufallsvariable Y := X +X +X reden, wobei angenommen sei, dass die drei VariablenX ,X ,X stochastisch unabhängig und identisch verteilt sind. Pseudozufallszahlen mit der Verteilung vonY lassen sich leicht generieren, etwa für exponentialverteilteXi durch > set.seed(687) > y < rexp(100) + rexp(100) + rexp(100) Was ist nun der Erwartungswert von Y, E(Y)? Wir können ihn bei dieser einfachen Funktion sofort angeben, denn der Erwartungswert ist additiv, der Erwartungswert für jedesXi ist , also ist E(Y) = . Nur ist die Simulation in diesem Fall nicht besonders effizient: > mean(y) [1] 2.828870 Was aber ist die Dichte vonY? Was ist die Quantilsfunktion? Oder: Was ist die Verteilungsfunktion? Wir können natürlich wieder mit hinreichend großem n Die Arithmetik von Zufallsvariablen ist sehr verschieden von der üblichen Arithmetik. Die Verteilung der SummeX +X +X ist i.d.R. nicht die Verteilung von X . . Funktionen von Zufallsvariablen simulieren und anschließend die empirische Dichte, Quantilsfunktion und Verteilungsfunktion benutzen. Das ist bei etwas komplizierteren Funktionen schnell sehr ineffizient. Zudem verhalten sich die Ergebnisse dieser Schätzer nicht so wie die d, q und p-Funktionen für bekannte Verteilungen. Oder man versucht, die Dichten analytisch zu finden. Das ist in diesem Beispiel zwar sehr einfach, im allgemeinen aber sehr mühsam. Das Paket distr erlaubt das symbolische Rechnen mit Funktionen von Zufallsvariablen. Insbesondere werden die d,p,q-Methoden für Funktionen von unabhängigen Zufallsvariablen bereitgestellt. Für unser Beispiel können wir damit schreiben: > library(distr) > X1 < Exp(1) > X1 Distribution Object of Class: Exp rate: 1 X1 ist ein Objekt, das eine exponentialverteilte Zufallsvariable mit Parameter (Rate) repräsentiert. Zufallsvariable können für alle im Basispaket stats vorhandenen Verteilungen erzeugt werden. Ihre Namen entsprechen den Namen imBasispaket ohne das Präfix d,p,q,r und beginnenmit einemGroßbuchstaben. Die Namen der Parameter sind in beiden Paketen bis auf wenige Ausnahmen gleich, ebenso wie deren Voreinstellungen. Mit diesen Objekten kann man nun durch Funktionen neue Objekte erzeugen, die wiederum entsprechende Zufallsvariable repräsentieren: > X1 < Exp(1) > X2 < Exp(1) > X3 < Exp(1) > Y < X1+X2+X3 > Y Distribution Object of Class: Gammad shape: 3 scale: 1 In diesem einfachen Fall ist sogar die exakte Verteilung identifiziert worden: Die Verteilung der Summe von drei exponentialverteilten Zufallsvariablen ist eine Gammaverteilung mit Formparameter . Versuchen wir ein komplizierteres Problem:Was sind Dichte, Verteilung und Quantilsfunktion der FunktionX exp(Y), wobeiX exponentialverteilt mit Parameter rate=2 ist undY unabhängig vonX normalverteilt mit Parametern mean=0.2, sd=0.5 ist? > X < Exp(rate=2) > Y < Norm(mean=0.2,sd=0.5) > Z < X∗exp(Y) > Z Distribution Object of Class: AbscontDistribution Grundlagen der Simulation Das ObjektZ repräsentiert eine absolut stetige Verteilung, eine Funktion der beiden ZufallsvariablenX undY. Für das Objekt Z existieren Methoden, die den d,p,q-Funktionen des Basispakets stats entsprechen bzw. sie approximieren. Für das Beispiel können wir schreiben: > d(Z)(1.3) # Dichte [1] 0.1799813 > p(Z)(1.3) # Verteilungsfunktion [1] 0.8510667 > q(Z)(0.5) # Median [1] 0.4099107 Es gibt zudem eine plot Methode für solche Objekte wie Z. Sie liefert einen Plot der Dichte, der Verteilungsfunktion und der Quantilsfunktion: > plot(Z) erzeugt Abbildung . . Abbildung 7.4: Dichte, Verteilungs- und Quantilsfunktion der Zufallsvariablen Z = X exp(Y), wennX exponentialverteilt mit Rate undY normalverteilt mit Erwartungswert . und Standardabweichung . ist. Das Paket distrEx erweitert einige der Möglichkeiten des Pakets distr. Insbesondere gibt es Funktionen für die Berechnung von Erwartungswerten der Zufallsvariablen sowie von Erwartungswerten von Funktionen der Zufallsvariablen. Die Notation folgt der üblichen Notation der Statistiker. Für die Verteilung der Variablen Z erhält man etwa: Die Grundlage der Methode ist eine Simulation der Funktionen mit den Methoden des Basispakets. Dichten, Quantile und Verteilungen werden dann automatisch mit den entsprechenden beschreibenden Methoden des letzten Kapitels berechnet. Allerdings werden Informationen über den Träger der Verteilung (die Menge der Zahlen positiver Wahrscheinlichkeit oder positiver Dichte) ebenso wie über den diskreten/absolut stetigen Charakter der Verteilung berücksichtigt. Zudem benutzt der + Operator, angewandt auf Objekte der Klasse distr, wenn möglich die schnelle Fourier-Transformation. Das Ergebnis ist in den meisten Fällen viel genauer als die direkte, nicht für den Anwendungsfall optimierte Simulation. . Übungsaufgaben > library(distrEx) > erw < E(Z) > erw [1] 0.6900549 > erw2 < E(Z, function(z) z^2) > erw2 [1] 1.201845 . Übungsaufgaben 1) Untersuchen Sie die Folge xi := ai ai , wobei a irrational ist. Zeichnen Sie insbesondere die aufeinanderfolgenden Werte xi, xi+ . Berechnen Sie die Korrelation zwischen x i und x i+ für i = , . . . , . 2) Untersuchen Sie die Folge xi := i π i π . Zeichnen Sie die Paare xi, xi+ . Berechnen Sie die zweiten Differenzen xi+ xi+ + xi und zeichnen Sie die ersten Werte. Was bedeutet das Ergebnis für die dreidimensionale Gleichverteilung der Folge? 3) Erzeugen Sie einen Vektor x, der Realisationen einer standardnormalverteilten Zufallsvariablen enthält. Wählen Sie einmal set.seed(123), dann set.seed(3). Berechnen Sie die Korrelation zwischen beiden Vektoren. 4) Benutzen Sie die Transformationsmethode aus Abschnitt . . , um Pseudozufallszahlen mit der Verteilung F(x) = x/( . + x) zu erzeugen. Berechnen Sie also die Quantilsfunktion F (p) und erzeugen Sie Pseudozufallszahlen durch F (Ui), wobeiUi, i = , . . . , gleichverteilte Pseudozufallszahlen sind. Vorsicht: Bei der Erzeugung der Pseudozufallsvariablen durch die Inverse der Verteilungsfunktion F (u) = . u/( u) muss im Zähler und im Nenner die gleiche Pseudozufallsvariable verwandt werden! Der Erwartungswert dieser Verteilung ist . Können Sie ein einfaches Argument finden, warum das so ist? Was passiert, wenn Sie die Anzahl der Pseudozufallsvariablen erhöhen? 5) Berechnen Sie Pr(X > ) = π( + x ) dx Die Variable hat also eine Cauchy-Verteilung. Benutzen Sie sowohl die pcauchy() Funktion und eineMonte-Carlo-Simulationmit der Funktion rcauchy(). Beachten Sie dabei, dass Pr(X > ) /n ni= 1(Xi > ) für n unabhängige Cauchy-VariableXi gilt. Grundlagen der Simulation 6) Die obige Cauchy-Verteilung ist symmetrisch um . Man kann das Integral daher als Pr(X > ) = π( + x ) dx schreiben. Das letzte Integral kann man wiederum als Erwartungswert einer Gleichverteilung auf dem Intervall [ , ] sehen, nämlich als E π( + U ) Berechnen Sie das Integral nochmal mit gleichverteilten Pseudozufallszahlen. Gibt es Unterschiede der Genauigkeit im Vergleich zu den Ergebnissen von pcauchy()? 7) Benutzt man das Paket distr, dann ist die Summe von unabhängigen und gleichverteilten Variablen etwa durch U < Unif() U2 < U+U U4 < U2+U2 U8 < U4+U4 U12 < U8+U4 6 gegeben. Die letzte Zeile normiert den Erwartungswert der Zufallsvariablen U12 auf . Diese Situation ist eines der am häufigsten benutzten Beispiele für den zentralen Grenzwertsatz: Die Summe unabhängiger, identisch verteilter Zufallsvariabler ist annähernd normalverteilt, zumindest wenn deren Varianz beschränkt ist. Vergleichen Sie die Dichte des Objektes U12 mit denWerten der Funktion dnorm() graphisch. Was ist die beste Wahl der Standardabweichung für die approximierende Normalverteilung? 8) Der folgende Code soll (inkorrekt) eine Gleichverteilung auf dem Einheitskreis simulieren: > set.seed(122) > r < runif(1000) > theta < runif(100,max=2∗pi) > plot(r∗cos(theta),r∗sin(theta)) > s < seq(0,2∗pi,by=0.05) ## Noch den Kreis > lines(cos(s),sin(s)) ## einzeichnen Warum funktioniert das nicht? Hinweis:Die Fläche eines Kreises mit Radius r ist πr . Also hat der Kreis mit Radius / eine Fläche von π/ . Der Kreisring (der Einheitskreis ohne den inneren Kreis mit Radius / ) hat also die Fläche π π/ = π/ , das dreifache der Fläche des inneren Kreises. Beide haben unter der obigen Simulation aber gleiche Wahrscheinlichkeit (warum?). Können Sie eine funktionierende Variante . Übungsaufgaben vorschlagen?Hinweis:Welche Funktion von r kannman nehmen, so dass die Wahrscheinlichkeit von Kreisen proportional zu ihrem Flächeninhalt wird? 9) Simulieren Sie eine Gleichverteilung auf der Kreisfläche, indem Sie zunächst je auf [- , ] gleichverteilte Pseudozufallszahlen X und Y erzeugen. Löschen Sie dann alle Paare, für die X^2 + Y^2> 1 ist. Warum ergibt sich eine Gleichverteilung? Benutzen Sie auch eine Simulation, um den Anteil an erzeugten Paaren zu schätzen, der gelöscht wird. 10) Aus einer Liste mit Elementen soll jedes Element mit Wahrscheinlichkeit p = . ausgewählt werden. Probieren Sie die folgenden drei Möglichkeiten und messen Sie den jeweiligen Zeitaufwand mit dem Befehl system.time() a) Erzeugen Sie n Bernoulli-Variable mit rbinom(). b) Die geometrische Verteilung ist die Anzahl von Versuchen, bis zum ersten Mal in einer Folge von unabhängigen Bernoulli-Variablen mit Wahrscheinlichkeit p eine auftaucht. Die Verteilung ist Pr(Y = k) = p( p)k. Man kann also die Indizes der Elemente, die ausgewählt werden, durch eine Summe von geometrisch verteilten Pseudozufallsvariablen erhalten. Der Befehl zur Erzeugung geometrisch verteilter Pseudozufallszahlen ist rgeom(). In dieser Version muss man im Durchschnitt nur np Zufallsvariablen generieren, statt der n in der vorigen Version. Allerdings muss man entweder sequentiell vorgehen oder von vornherein eine entsprechend große Anzahl geometrischer Variabler generieren. c) Man bestimmt zunächst die Anzahl der auszuwählenden Elemente durch den Wert m einer binomialverteilten Pseudozufallsvariablen mit Parametern n und p. Dann benutzt man die Funktion sample(n,m). Warum funktioniert diese Variante?

Chapter Preview

References

Zusammenfassung

Vorteile

- Einführung in die statistische Analyse mit R für Wirtschafts- und Sozialwissenschaftler

- Inklusive hilfreicher Tipps wie "Ansprechende Grafiken mit R gestalten"

Zum Thema

R ist ein Statistikprogramm, das kostenlos über das Internet verbreitet wird und dessen Source Codes frei zugänglich sind.

Aufgrund dieses kostenlosen Angebots gehen immer mehr Dozenten dazu über, neben SPSS auch R zu lehren bzw. SPSS durch R zu ersetzen.

In R steht dem Nutzer die gesamte Bandbreite statistischer Verfahren zur Verfügung. Durch die eigenständige Programmierumgebung ist die Software sehr flexibel und erlaubt notwendige Modifikationen und Erweiterungen verfügbarer Prozeduren.

Zum Werk

Dieses Buch führt leicht verständlich in die statistische Analyse mit R ein. Anhand von Beispielen wird die Umsetzung der wichtigsten Methoden der Statistik, wie sie üblicherweise in den Grundkursen gelehrt werden, mit R vorgestellt.

Das Buch verfolgt entsprechend zwei Ziele:

1. Vorstellung der statistischen Methoden,

2. Benutzung des Werkzeuges R zur Analyse von Daten.

Inhalt

- Grundlagen von R

- Datenbehandlung und graphische Darstellungen mit R

- Datenbeschreibungen (deskriptive Statistik)

- Wahrscheinlichkeitsverteilungen

- Regressionsanalysen

- Optimierungsverfahren

- Simulationen mit R

Neben vielen neuen, wirtschaftsorientierten Beispielen wird nun auch in die Paneldatenanalyse und Stichprobentheorie eingeführt.

Zu den Autoren

Dr. Andreas Behr ist wissenschaftlicher Mitarbeiter am Institut für Statistik und Ökonometrie der Universität Münster.

Dr. Ulrich Pötter ist wissenschaftlicher Mitarbeiter am Institut für Statistik der Universität Bochum.

Zielgruppe

Für Studierende und Dozenten der Wirtschaftswissenschaften im Bachelor an Universitäten und Fachhochschulen.