S M Simulierte Daten geben offenbar keine Auskunft über die Realität. Aber sie ermöglichen es, Konsequenzen theoretischer Modelle mit realen Daten zu vergleichen. Und man kann Eigenschaften statistischer Algorithmen in unterschiedlichen, aber wohldefinierten Situationen ausprobieren. Diese beiden Aspekte untersuchen wir mit den Möglichkeiten, die R bereitstellt. Zunächst wird das Standardmodell statistischer Argumentation eingeführt und entsprechende R Befehle benutzt, um diese Argumente nachzuvollziehen. Dann werden einige einfache wahrscheinlichkeitstheoretische Modelle abhängiger Variabler eingeführt, deren Implementation in R vorgestellt und ihr potentieller Nutzen für statistische Anwendungen illustriert. 8.1 Das Standardmodell der Statistik 8.1.1 Mittelwert oder Median? 8.1.2 Warum n ? 8.1.3 Kann man Erwartungswerte schätzen? 8.1.4 Bootstrap und Bagging 8.2 Markow-Ketten 8.3 Übungsaufgaben . Das Standardmodell der Statistik Der Ausgangspunkt für die beschreibenden Methoden in den Kapiteln und war die Annahme, man könne für eine vorgegebene Liste von identifizierten Personen (oder Organisationen, Institutionen, Ereignissen etc.) jeder Person (i.A. jeder identifizierten Organisation etc.) einen Wert x aus einer vorab bestimmten Menge eindeutig zuordnen. Daraus ergab sich die Aufgabe, die Verteilung dieser Werte sowie deren Eigenschaften angemessen zu repräsentieren. Bei dieser Fragestellung der beschreibenden Statistik gibt es keinen Platz für Wahrscheinlichkeitstheorie oder gar Simulationsverfahren. Die verschiedenen beschreibenden Verfahren müssen im jeweiligen Verwendungskontext beurteilt werden. Ein Verfahren ist besser als ein anderes Verfahren, wenn es Stochastische Modelle für die gegebenen Daten und den entsprechenden Kontext die „Information“ für potenzielle Nutzer „besser“ (präziser, kompakter, einfacher, übersichtlicher, verständlicher, schöner,. . . ) darstellt. Solche allgemeinen, kontext- und nutzerabhängigen und wenig präzisen Kriterien sind aber kaum geeignet, allgemeine Empfehlungen für bestimmte Verfahren zu entwickeln. Ein möglicher einheitlicher Rahmen für die allgemeine Beurteilung statistischer Verfahren wurde erst seit Anfang des letzten Jahrhunderts entwickelt. Er verbindet Wahrscheinlichkeitstheorie und statistische Praxis durch ein einfaches wahrscheinlichkeitstheoretisches Modell, das Aussagen über die „Güte“ von statistischen Verfahren ohne Bezug auf gegebene Daten und deren Kontexte zulässt. Die Idee ist einfach, aber radikal. Der Bezug zu realen Einheiten u und deren zugewiesenen Werten X(u) wird ersetzt durch die Betrachtung von unabhängigen und identisch verteilten Zufallsvariablen. Die Menge wird nur noch zur Indizierung eines Vektors von unabhängigen Kopien von ZufallsvariablenX : Ω benutzt. Statistische Verfahren werden dann als Funktionen dieser Zufallsvariablen interpretiert. Damit kann man Methoden der Wahrscheinlichkeitstheorie verwenden, um für große Klassen von Verteilungen die Optimalität von bestimmten statistischen Verfahren zu beweisen. Die Optimalität hängt dann gar nicht mehr von den spezifischen Aspekten vorliegender Daten ab, sie gilt für eine sehr große, wenn auch beschreibbare Klasse von Daten. Die Simulationsmöglichkeiten des letzten Kapitels machen auch klar, wie groß diese Klasse ist: Sie umfasst alle Datensätze, die sich durch Simulation mit vorgelegten Verteilungen erstellen lassen. Heute ist die Einbettung statistischer Verfahren in ein wahrscheinlichkeitstheoretischesModell so verbreitet, dassmanchmal sogar der Eindruck vermittelt wird, Daten müssten auch wie in dieser Interpretation erzeugt sein, um berechtigt statistische Verfahren benutzen zu können. Man übersieht dabei, dass Mathematik nichts über reale Daten und ihre Bezüge wissen kann (es auch gar nicht will), folglich auch nicht vorschreiben kann, wie man mit Daten umgeht. In der Tat ist der Nutzen des Modells und damit der Einführung allgemeiner Kriterien für die Beurteilung statistischer Verfahren viel indirekter. Es muss gar nicht vorausgesetzt werden, Daten seien als Realisation unabhängiger und identisch verteilter Zufallsvariabler aus einer Klasse von Verteilungen entstanden. Denn wenn ein Verfahren für alle unter dem mathematischen Modell erzeugbaren Daten schlechter ist als ein anderes, dann bedarf es wohl kontextbezogener Gründe, dennoch das schlechtere Verfahren für die vorliegenden Daten zu wählen. . . Mittelwert oder Median? Ein klassisches Beispiel ist die Wahl zwischen Mittelwert und Median als Maßzahl für die Lage eines Datensatzes. Wir übersetzen das Problem zunächst in ein wahrscheinlichkeitstheoretisches Modell und dann in ein . Das Standardmodell der Statistik entsprechendes R Programm: Wir betrachten n (unabhängige und identisch verteilte) Zufallsvariable mit vorgegebener Verteilung, etwaX , . . . ,Xn und definieren Mittelwert und Median wie in Kapitel : M(X , . . . ,Xn) := n n i= Xi , med(X , . . . ,Xn) := F̂ - X ,...,Xn( / ) Allerdings entsprechen die Ausgangsdaten nun nicht einer Liste von Werten, sondern sie werden als Realisationen von Zufallsvariablen aufgefasst. Mittelwert und Median sind daher als Funktionen von Zufallsvariablen ebenfalls zufällige Größen. In R wählen wir als Verteilung eine Standardnormalverteilung und setzen n = . Dann ergeben sich für zwei Realisationen etwa: > set.seed(132) > n < 10 > x < rnorm(n) > mean(x) [1] 0.08055975 > median(x) [1] 0.1653128 > ### Neue Ziehung: > x < rnorm(n) > mean(x) [1] 0.03467084 > median(x) [1] 0.1580355 Nun sind die Ausgangswerte für die Simulation bekannt: Der Erwartungswert und der Median der Standardnormalverteilung sind beide . Die Funktionen Mittelwert und Median, M(X , . . . ,X ) und med(X , . . . ,X ), werden nun als Schätzfunktionen des zugrundeliegendenWertes aufgefasst undman kann fragen, wie nah diese Schätzfunktionen diesemWert kommen. Am einfachsten lässt sich der erwartete quadratische Abstand handhaben. Gesucht ist also der Wert von E (M(X , . . . ,X ) ) bzw. E (med(X , . . . ,X ) ) wobei sich der Erwartungswert auf die Verteilung der Zufallsvariablen (X , . . . , X ) bezieht. Zwar kann man diese Größen auch theoretisch berechnen, hier aber benutzen wir eine Simulation, um die Erwartungswerte näherungsweise zu berechnen. Wir benutzen Simulationen und tragen die Ergebnisse in die Matrix erg ein: > set.seed(132) > n < 10 > nsim < 100000 Stochastische Modelle > erg < matrix(0,nrow=nsim,ncol=2) > for(i in 1:nsim){ > x < rnorm(n) > erg[i,1] < mean(x) > erg[i,2] < median(x) } Zur Auswertung berechnen wir den Durchschnitt der (erg[i,] 0)^2 über alle Simulationen. Dazu benutzen wir den Befehl colMeans() der die Mittelwerte getrennt für die Spalten der Matrix erg berechnet. > colMeans(erg^2) [1] 0.09882695 0.13703214 Der Median erweist sich gegenüber demMittelwert als deutlich schlechter, zumindest wenn als Kriterium die erwartete quadratische Abweichung gewählt wird. Die erwarteten quadrierten Abstände des Medians sind fast % größer als die des Mittelwerts. Man sagt auch, der Mittelwert sei effizienter als der Median. Dieses Effizienzkriterium lässt sich auch direkt in ein Erfordernis an benötigte Fallzahlen übersetzen: Der erwartete quadratischen Abstand des Mittelwerts bei Beobachtungen ist . . Um diesen Wert mit dem Median als Schätzfunktion zu erreichen, benötigt man % mehr Beobachtungen, also etwa . Man kann auch zeigen, dass dieses Verhältnis nicht von der Fallzahl n abhängt. . . Warum n ? Wir können nun auch ein Argument nachvollziehen, das zu der Verwendung des Faktors /(n ) an Stelle von /n führt. Man argumentiert, dass die Schätzfunktion var(X , . . . ,Xn) := /n n i (Xi M(X , . . . ,Xn)) im Durchschnitt zu klein ist. Denn der Mittelwert M(.) minimiert gerade die Quadratsumme var(.), der Wert ist also kleiner als der Wert, der sich ergeben würde, wenn man an Stelle des Mittelwerts den Erwartungswert der zugrundeliegenden Wahrscheinlichkeitsverteilung verwenden würde. In der Tat ist E var(X , . . . ,Xn) = E n n i= Xi M(X , . . . ,Xn) = V(X ) wobei V(X) := E(X ) (E(X)) die Varianz der zugrundeliegenden Verteilung ist. In anderen Worten: Für alle möglichen Werte der Varianz der zugrundeliegenden Verteilung ist der Erwartungswert der Schätzfunktion var(.) (mit dem Faktor /(n )) gleich der Varianz. Man sagt auch, die Schätzfunktion sei erwartungstreu. Wir simulieren das für n = : > set.seed(456) Der entsprechende Befehle für Zeilen einer Matrix ist rowMeans(). Diese Befehle sind etwas effizienter als die Äquivalente apply(erg,1,mean) bzw. apply(erg,2,mean). . Das Standardmodell der Statistik > n < 5 > nsim=10000 > erg < matrix(0,ncol=2,nrow=nsim) > for(i in 1:nsim){ + x < rnorm(n) ##Standardnormalverteilung + erg[i,1] < var(x) + erg[i,2] < mean((x mean(x))^2) } > colMeans(erg) [1] 1.0087131 0.8069705 Der Durchschnitt der Werte von var(.) (als Simulationsapproximation des Erwartungswerts) liegt nah bei . Dagegen unterschätzt die Variante mit dem Faktor /n die Varianz um / . Die Eigenschaft der Erwartungstreue lässt sich offenbar nur im wahrscheinlichkeitstheoretischen Modell der Statistik formulieren, denn sie setzt eine zugrundeliegende Verteilung und deren Kennwerte (hier: die Varianz) immer schon voraus. Aber kann man daraus ableiten, man solle auch immer die erwartungstreue Variante var(.) benutzen? Was passiert etwa, wenn man sich für die Standardabweichung und nicht direkt für die Varianz interessiert? In R ist die entsprechende Schätzfunktion sd() als Wurzel der Funktion var() implementiert. Simulieren wir ihr Verhalten: > set.seed(456) > n < 5 > nsim=10000 > erg < matrix(0,ncol=1,nrow=nsim) > for(i in 1:nsim){ + x < rnorm(n) + erg[i,1] < sd(x) } > colMeans(erg) [1] 0.9435652 Der Durchschnitt ist nun kleiner als die zugrundeliegende Standardabweichung, diese Schätzfunktion ist also nicht mehr erwartungstreu. . . Kann man Erwartungswerte schätzen? Kehren wir noch einmal zumMittelwert und dessen Eigenschaften als Schätzfunktion zurück. Der Mittelwert hat zwar eine kleinere erwartete quadrierte Abweichung vom Erwartungswert als der Median in unserem Beispiel. Aber er Eine erwartungstreue Version müsste den unschönen Faktor Γ((n )/ )/( Γ(n/ )) verwenden, vgl. Erich Lehmann/George Casella:Theory of Point Estimation, Springer , S. . Diese Formel ist wohl noch nie in der angewandten Literatur benutzt worden. Das liegt nicht allein an der komplizierten Form, schließlich könnte man die Formel in einem Befehl eines Statistikprogramms wie R berechnen lassen. Das Problem ist vielmehr, dass der einfache Bezug zur beschreibenden Statistik verloren geht. Stochastische Modelle besitzt auch eine beunruhigende Eigenschaft: Er ist sehr sensitiv gegenüber Ausreißern. Einige wenige ungewöhnliche Werte können die Ergebnisse beliebig beeinflussen. Wir können das mit Hilfe des wahrscheinlichkeitstheoretischen Modells nun genauer formulieren. Dazu verändern wir die zugrundeliegende Verteilung um ein weniges, etwa indem wir mit sehr kleinerWahrscheinlichkeit einen sehr großen Wert setzen, sonst aber die unveränderte Verteilung zu Grunde legen.Hat alsoX die Verteilung F(.), dann setzen wir X := ( Y)X +Yx für eine vonX unabhängige Bernoulli-VariableY mit E(Y) = p. Der Erwartungswert vonX ist dannE(X ) = ( p)E(X)+px . Je nachWahl von p und x können sich die Erwartungswerte vonX undX beliebig unterscheiden, aber die Verteilungsfunktionen sind sich bei sehr kleinem p sehr ähnlich. Anders ausgedrückt: Der Erwartungswert verändert sich nicht stetig mit der zugrundeliegenden Verteilungsfunktion. Wir simulieren dieses Phänomen und setzen p = / und x = . Ist die Verteilung vonX wieder die Standardnormalverteilung, dann ist der Erwartungswert E(X ) gerade . > set.seed(123) > n < 10 > p < 0.00001 > x0 < 100000 > nsim < 10000 > erg < matrix(0,ncol=2,nrow=nsim) > for (i in 1:nsim){ + x < rnorm(n) + xs < ifelse(runif(n) colMeans(erg) [1] 2.0030888343 0.0008036376 > mean((erg[,1] 1)^2) [1] 19997.42 > mean(erg[,2]^2) [1] 0.1376943 Der Mittelwert der xs über Simulationen ist . Zudem liegen die meisten Ergebnisse bei , in einigen wenigen Fällen aber in der Nähe von . Weder diese beiden häufigsten möglichen Ergebnisse noch der Durchschnitt über die Simulationen sind auch nur in der Nähe des tatsächlichen Erwartungswerts. Der mittlere quadratische Abstand beträgt , so dass die Schätzfunktion sehr instabil ist. Weder eine Erhöhung von n noch die Wahl einer anderen Verteilung können dieses Problem beseitigen, denn man kann immer ein entsprechend kleineres p und gleichzeitig ein größeres x wählen. Kann man . Das Standardmodell der Statistik also nicht von vornhinein solche Verteilungen ausschließen, dann ist die Schätzfunktion Mittelwert instabil. Es kann daher in Abhängigkeit von der Klasse der möglichen Verteilungen beliebig schwer sein, Erwartungswerte zu schätzen. Der Median hat dieses Problem in weit geringeremMaße, er hängt nicht von den Rändern der Verteilung ab. Allerdings fallen nun Median und Erwartungswert der zugrundeliegenden Verteilung nicht mehr zusammen. Der Median der Verteilung vonX ist etwas größer als , der Erwartungswert ist . Die Schätzfunktion Median aber ist auch in diesem Beispiel robust: Der Mittelwert über alle Simulationen ist , und der mittlere quadratische Abstand beträgt nur . . . . Mittelwert oder Median: Bootstrap und Bagging Kann man nicht die guten Eigenschaften von Median und Mittelwert kombinieren und eine Schätzfunktion konstruieren, die effizienter als der Median und robuster als der Mittelwert ist? Dafür gibt es in der Tat viele Vorschläge. Eine einfache Idee der Kombination besteht darin, wiederholt den Median von allen möglichen Stichproben aus den vorgelegten Daten zu berechnen und aus diesen Werten den Mittelwert zu berechnen. Das Verfahren, aus vorgegebenen Daten Stichproben zu ziehen und für diese Stichproben wiederholt Schätzfunktionen zu berechnen, nennt man Bootstrap. Das Verfahren wird i.d.R. benutzt, um die Variabilität von Schätzfunktionen zu berechnen bzw. um Konfidenzintervalle und Tests zu konstruieren. Wir wollen es benutzen, um Schätzfunktionen zu kombinieren. In diesem Zusammenhang nennt man die Kombination auch Bagging (ein Akronym für „Bootstrap Aggregation“). Um das Verfahren durchzuführen, schreiben wir uns eine Funktion, in der zunächst nboot Stichproben (mit Zurücklegen) aus den Daten x gezogen werden. Anschließend wird für jede Stichprobe der Median berechnet und dann der Mittelwert aus allen Medianen gebildet: > ###Median Bagging > baggedmed < function(x,nboot){ + erg < vector("numeric",nboot) + n < length(x) Darauf haben Bahadur und Savage schon hingewiesen. Sie zeigen, dass es keine vernünftigen Tests, Konfidenzintervalle oder Schätzfunktionen für den Erwartungswert geben kann, wenn man nur eine hinreichend große Klasse von Verteilungen zulässt (R.R. Bahadur, Leonard J. Savage:The nonexistence of certain statistical procedures in nonparametric problems. Annals of Mathematical Statistics , , S. – ). Der Stichprobenumfang kann natürlich variiert werden. I.d.R. verwendetman aber den gleichen Umfang wie die Anzahl der vorgelegten Daten. Was bei der Variation der Stichprobenumfänge passiert und wie sich der Median von wiederholt berechneten Mittelwerten verhält, diskutiert Jose R. Berrendero:The bagged median and the bragged mean.The American Statistician, , , S. – . Stochastische Modelle + for(i in 1:nboot){ + erg[i] < median(sample(x,size=n,replace=T))} + mean(erg) } Wir probieren das Verfahren zunächst bei der Standardnormalverteilung aus. Wir müssen noch überlegen, wie oft Stichproben gezogen werden sollen. Da jeweils der Median berechnet wird, können nur wenige verschiedene Werte auftreten. Eine relativ geringe Anzahl von Stichproben sollte daher schon ausreichen. Wir wählen nboot < 20: > set.seed(132) > nboot < 20 #Bootstrap Stichproben > nsim < 10000 > n < 10 > ### > erg < vector("numeric",nsim) > for(i in 1:nsim){ + x < rnorm(n) + erg[i] < baggedmed(x,nboot) } > mean(erg) [1] 0.001678400 > mean(erg^2) [1] 0.1212633 Der Mittelwert über die Simulationen ist praktisch , der erwartete quadrierte Abstand ist . , also kleiner als die . des Medians, aber größer als die . des Mittelwerts. Ist diese Kombination von Median und Mittelwert auch noch robust? Wir simulieren das Beispiel aus dem letzten Abschnitt mit dem „bagged“ Median: > p < 0.00001 > x0 < 100000 > set.seed(132) > nboot < 20 #Anzahl Bootstrap Ziehungen > nsim < 10000 > n < 10 > erg < vector("numeric",nsim) > for(i in 1:nsim){ + x < rnorm(n) + y < ifelse(runif(n) mean(erg) [1] 0.002405802 > mean(erg^2) Bei ungerader Fallzahl und voneinander verschiedenen Werten können nur die vorliegenden n Werte auftreten, bei gerader Fallzahl auch noch die Mittelwerte von je zwei verschiedenen Werten. . Markow-Ketten [1] 0.1239421 Der Mittelwert über die Simulationen ist praktisch , entspricht also dem Median der zugrundeliegendenVerteilung. Der erwartete quadratische Abstand (zu ) ist wieder . , also immer noch besser als der einfache Median. . Markow-Ketten Bisher haben wir nur unabhängige und identisch verteilte Pseudozufallszahlen betrachtet. In vielen Fällen ist es aber wichtig, die Entwicklung von Größen über die Zeit oder Abhängigkeiten zwischen Ereignissen oder zwischen räumlich benachbarten Größen zu simulieren. In diesem Abschnitt besprechen wir einige einfache Möglichkeiten, solche Simulationen mit R durchzuführen. Die wohl einfachste Form der Abhängigkeit ergibt sich, wenn die weitere Entwicklung des Prozesses nicht von der gesamten Vorgeschichte abhängt, sondern nur vom zuletzt erreichten Wert. Repräsentiert := , , , . . . eine Zeitachse und t einen Zeitpunkt, dann soll die Folge von diskreten ZufallsvariablenX ,X ,X , . . . eineMarkow-Kette heißen, wenn Pr(Xt = xt Xt = xt , . . .X = x ) = Pr(Xt = xt Xt = xt ) gilt. Ist Pr(Xt = xt Xt = xt ) unabhängig von t, dann spricht man von einer homogenenMarkow-Kette. Einige Beispiele sollen das Konzept verdeutlichen. Ein Zustandsraum mit nur zwei Zuständen, etwa = , , ist die einfachste Möglichkeit. Man kann z.B. an die Zustände „arbeitslos“ und „beschäftigt“ denken. Eine homogene Markow-Kette mit diesem Zustandsraum ist vollständig durch den Startwert X = x (oder die Startverteilung Pr(X = x )) und die vier bedingten Wahrscheinlichkeiten Pr(Xt = Xt = ), Pr(Xt = Xt = ), Pr(Xt = Xt = ) und Pr(Xt = Xt = ) beschrieben. Arrangiert man diese Wahrscheinlichkeiten in einer Matrix Pt := Pr(Xt = Xt = ) Pr(Xt = Xt = ) Pr(Xt = Xt = ) Pr(Xt = Xt = ) dann lassen sich die marginalen Wahrscheinlichkeiten einer homogenen Markow-KetteXt , bei der Pt ja für alle Zeitpunkte gleich einer Matrix P ist, als Pr(Xt = .) = (Pr(X = ), Pr(X = )) P P P = (Pr(X = ), Pr(X = )) Pt Wir werden gleich sehen, dass sich das Konzept für beliebige Wertebereiche formulieren lässt. Die Beschränkung auf diskrete endliche Zufallsvariable vereinfacht nur die Formulierung der Markow-Bedingung. Zudem benutzt man im diskreten Fall andere Simulationstechniken als im stetigen Fall. Stochastische Modelle schreiben, wobei „*“ für das Matrixprodukt steht. Nun sind wir an der Simulation solcher Prozesse interessiert, nicht nur an der Verteilung an bestimmten Zeitpunkten. Dazu können wir die Indexschreibweise von R nutzen. Ist etwa P = . . . . undX = , dann lässt sich eine Realisation der ersten Werte dieser Markow- Kette durch > set.seed(823) > P < matrix(c(0.4,0.6,0.2,0.8),nrow=2,ncol=2,byrow=T) > Schritte < 20 > erg < vector("numeric",Schritte+1) > erg[1] < 1 ##Startwert > for(i in 2:(Schritte+1)){ + erg[i] < sample(1:2,size=1,prob=P[erg[i 1],])} > erg [1] 1 2 2 2 2 2 1 1 2 2 2 1 1 1 1 1 2 2 2 2 1 realisieren. Benutzt wird nur der Befehl sample() in Verbindung mit der vorgegebenen Übergangsmatrix P und der Indexmöglichkeit von R. Markow-Ketten konvergieren unter bestimmten Regularitätsannahmen zu einer stationären Verteilung: Unabhängig vom Startwert x ist die Verteilung vonXt für große t konstant mit einer Verteilung, die nur von den Eigenschaften von P abhängt. Im obigen Beispiel gilt für die Verteilung Pr(X = ) = . und Pr(X = ) = . gerade ( . , . ) P = ( . , . ) es ist also die eindeutige stationäre Verteilung der Markow-Kette. Wir überprüfen das in R: > b < c(0.25,0.75) > b%∗%P==b [,1] [,2] [1,] TRUE TRUE Wir können unsere Simulation etwas länger laufen lassen und die marginale Verteilung nach Schritten betrachten. Sie sollte schon nah an der stationären Verteilung liegen: > set.seed(823) > Schritte < 1019 > erg < vector("numeric",Schritte+1) Eine gute Einführung in dieTheorie der Markow-Ketten und eine ausführliche Diskussion von Regularitätsannahmen findet sich bei James R. Norris: Markov Chains, Cambridge University Press, . . Markow-Ketten > erg[1] < 1 ##Startwert > for(i in 2:(Schritte+1)){ + erg[i] < sample(1:2,size=1,prob=P[erg[i 1],])} > prop.table(table((erg[21:(Schritte+1)]))) 1 2 0.257 0.743 In vielen Fällen lassen sich Markow-Ketten als iterierte zufällige Funktionen generieren. Das ist auch die Form, in der man Markow-Ketten mit beliebigen Zustandsräumen simulieren kann. Denn dann ist die Anzahl der möglichen Übergänge i.d.R. viel zu groß, um noch effizient durch eine Matrix dargestellt werden zu können. Man erhält z.B. eine Irrfahrtauf den ganzen Zahlen . . . , , , , . . ., wenn man etwaX = setzt und dannXt := Xt berechnet, wobei + bzw. mit gleichen Wahrscheinlichkeiten gewählt werden. Etwas allgemeiner kann manXt := aXt wählen und zudem an Stelle der Gleichverteilung auf , eine andere Verteilung verwenden. > set.seed(213) > a < 0.5 #Faktor > p < 0.5 #Wahrscheinlichkeit für +1 > schritte < 250 > start < 0.3 > zuf < 2∗(runif(schritte) Markow < function(x,y) a∗x+y > erg < Reduce(Markow,zuf,init=start,accumulate=T) Abbildung 8.1: Irrfahrt mit a = . . Der Befehl Reduce() hat als erstes Argument eine Funktion mit zwei Argumenten. Diese Funktion wird rekursiv auf die Elemente der Liste oder des Vektors angewandt, der als zweites Argument angegeben werden muss. Das dritte Argument ist der Startwert. Berechnet wird also erst Markow(start,zuf[1]), dann Markow(Markow(start,zuf[1]), zuf[2])usw.Mit derOption accumulate =T werden die Zwischenergebnisse ebenfalls ausgegeben. Die Abbildung . zeigt das Ergebnis. Was ist die stationäre Verteilung dieser Kette? Wir starten Markow- Ketten mit einer Normalverteilung mit Standardabweichung und berechnen einen Dichteschätzer nach Schritten: > schritte < 100 > nsim < 100000 Stochastische Modelle > start < 5∗rnorm(nsim) > zuf < as.data.frame(matrix(2∗(runif(schritte∗nsim) erg < Reduce(Markow,zuf,init=start,accumulate=T) > plot(density(erg[[schritte+1]])) Abbildung 8.2: Dichte der stationären Verteilung der Irrfahrt mit a = . . Es ergibt sich offenbar eine Gleichverteilung auf dem Intervall ( , ) (Abbildung . ). Das kann man in diesem einfachen Fall auch direkt sehen: IstXt gleichverteilt auf ( , ), dann istXt entwederXt / + oder Xt / . Der erste Ausdruck ist gleichverteilt auf ( , ), der zweite gleichverteilt auf ( , ), und eine von beiden Möglichkeiten wird mit Wahrscheinlichkeit / gewählt. Die Verteilung vonXt ist also wieder die Gleichverteilung auf ( , ). Die zufällige Transformation ändert die Gleichverteilung nicht. Folglich ist es die stationäre Verteilung der Markow-Kette. Markow-Ketten, die als iterierte zufällige Funktionen konstruiert werden, haben eine eigentümliche Eigenschaft, die man gut bei Simulationen ausnutzen kann. Während die Kette wie in Abbildung . den Zustandsraum ganz irregulär durchläuft, konvergiert die Kette sehr schnell, wenn man die zufälligen Funktionen „rückwärts“ anwendet. Sind f (.), f (.), . . . , fm(.) die zufälligen Funktionen, dann hatten wir die Markow-Kette durch x , f (x ), f (f (x )), . . . konstruiert. Beim Rückwärtsprozess wendet man die Funktionen in umgekehrter Reihenfolge an, also als x , f (x ), f (f (x )), f (f (f (x ))), . . .. In R kann man sich eine entsprechende Funktion definieren, die diesen Rückwärtsprozess simuliert. Dazu berechnen wir zunächst x , f (x ), f (x ), . . ., dann wenden wir f (.), f (.), . . . auf den zuvor berechneten Vektor x , f (x ), f (x ), . . . ohne die beiden ersten Elemente an und wiederholen das Verfahren bis zur gewünschten Schrittzahl: > rueckwaerts < function(start,zuf,f=Markow){ + n < length(zuf) Einzelheiten und weitere Beispiele finden sich bei Persi Diaconis, David Freedman: Iterated random functions. SIAM Review, , , S. – . Weitere graphische Darstellungen iterierter zufälliger Funktionen sind auf den Seiten – nach dem Aufsatz dargestellt. Eins dieser Beispiele wird in der Aufgabe ) zu diesem Kapitel aufgegriffen. In der Zeitreihenanalyse werden allgemeine Irrfahrten mit beliebigem a und beliebigen (aber zeitkonstanten) Verteilungen des additiven Terms autoregressive Prozesse genannt. Eigenschaften der stationären Verteilung sind selbst für diese einfachen Irrfahrten erstaunlich schwer zu bestimmen, wenn a = . bzw. a = ist. Im Aufsatz von Diaconis/Freedman finden sich einige Hinweise in den Abschnitten . und . . . Markow-Ketten + erg < vector("numeric",n+1) + ind < 1:n + erg[1] < start + erg[1+ind] < f(start,zuf)#f_1(start),f_2(start),... + for(i in 1:(n 1)){ + ind < ind[ length(ind)] + erg[1+i+ind] < f(erg[1+i+ind],zuf[ind])} + return(erg) } (a) (b) Abbildung 8.3: Rückwärtsprozess der Irrfahrt mit a = . . a) Die ersten Schritte mit den zufälligen Funktionen aus Abbildung . . b) verschiedene Startwerte mit gleichen zufälligen Funktionen. Abbildung 8.4: Dichte des Rückwärtsprozesses der Irrfahrt mit a = . nach Schritten. Wir simulieren nun den Rückwärtsprozess mit den gleichen Zufallszahlen wie zuvor und zeichnen die ersten Schritte: > set.seed(213) > schritte < 20 > start < 0.3 > zuf < 2∗(runif(schritte) plot(rueckwaerts(start,zuf), + type="o") Der Grenzwert wird in nur Schritten erreicht und hat den Wert . . Hängt dieser Grenzwert vom Startwert x ab? Wir versuchen verschiedene Startwerte, aber benutzen in allen Fällen die gleichen zufälligen Funktionen. Das Ergebnis ist in Abbildung . dargestellt: Die Rückwärtsprozesse konvergieren unabhängig vom Startwert zum gleichen Wert. Der Wert hängt aber von den zufälligen Funktionen f (.), f (.), . . . ab. Die Grenzverteilung ist Stochastische Modelle die stationäre Verteilung des Vorwärtsprozesses. Um das zu sehen, simulieren wir Rückwärtsprozesse mit jeweils Schritten und zeichnen einen Kerndichteschätzer. Das Ergebnis zeigt Abbildung . . > set.seed(213) > schritte < 20 > p < 0.5 > a < 0.5 > nsim < 10000 > zuf < matrix(2∗(runif(nsim∗schritte) erg < matrix(0,nrow=nsim,ncol=schritte+1)#Startwert 0 > for(i in 1:nsim) { + erg[i,1:(schritte+1)]< rueckwaerts(erg[i,1],zuf[,i])} > plot(density(erg[,schritte+1])) . Übungsaufgaben 1) Benutzen Sie eine Simulation, um die relative Effizienz des Mittelwerts gegenüber dem Median für eine t-Verteilung mit Freiheitsgraden und n = zu berechnen. Hinweis: Pseudozufallszahlen der t-Verteilung erhält man mit rt(n,df=5). Was passiert, wenn Sie df=1.5 wählen? 2) Überlegen Sie, warum man die Simulationsergebnisse in der vorherigen Aufgabe (ebenso wie für die Normalverteilung in Abschnitt . . ) nicht für weitere Parameterwerte von Median bzw. Erwartungswert der zugrundeliegenden Wahrscheinlichkeitsverteilung wiederholen muss. Hinweis: Wie verändern sich die Werte der Schätzfunktionen, wenn zu allen Werten der ZufallsvariablenX eine Konstante addiert wird? Wie verändern sich dann Erwartungswert und Median der zugrundeliegenden Wahrscheinlichkeitsverteilung? Und wie die Werte der Kriteriumsfunktion erwartete quadrierte Abweichung? 3) Ein Modell diffusiver Aggregation: Eine x Matrix A enthält zu Anfang die Werte , nur die letzte Zeile enthält . Dann werden nacheinander „Partikel“ gestartet. Sie beginnen in einer zufälligen Position in der ersten Zeile (benutzen Sie sample()). Sie bewegen sich bei jedem Schritt mit den Wahrscheinlichkeiten ( . , . , . , . ) nach oben, nach unten, nach rechts oder nach links (eine zweidimensionale Irrfahrt). Die Irrfahrt endet, wenn entweder das Partikel die Matrix verlässt Das wird in der „perfekten Simulation“ (oder: „coupling from the past“) ausgenutzt, vgl. James Propp, David Wilson: Coupling from the past: A user’s guide. Microsurveys in discrete probability, DIMACS Ser. Discrete Math. Theoret. Comput. Sci., , , S. – . Persi Diaconis:The Markov chain Monte Carlo Revolution, Bulletin of the American Mathematical Society, , , S. – , enthält weitere Beispiele. . Übungsaufgaben oder wenn es sich in der Nachbarschaft eines Matrixelements mit dem Wert befindet (genauer: wenn es sich genau eine Zeile über bzw. unter einem Element bzw. eine Spalte rechts oder links von einem Element mit demWert befindet). In letztem Fall wird die letzte Position des Partikels in der Matrix auf gesetzt. Starten Sie solange Partikel bis insgesamt Werte auf gesetzt wurden. Malen Sie das Ergebnis mit image(A,col=c("white","black"). Ist die Folge der Matrizen A$_{t}$ eine Markow-Kette? Ist der Prozess homogen? 4) Starten Sie einen Prozess an der Stelle x des Intervalls ( , ), etwa . . Wählen Sie zufällig eine Richtung: nach links (kleiner als x ) oder nach rechts (größer als x ) mit Wahrscheinlichkeit p = / . Wählen Sie dann je nach Wahl der Richtung einen Punkt y gleichverteilt aus dem Intervall ( , x ) oder aus (x , ). Wiederholt man das Vorgehen, erhält man eine Markow-Kette mit Zustandsraum ( , ). Schreiben Sie ein Programm, das Schritte der Markow-Kette simuliert. Hinweis: Ist U eine auf ( , ) gleichverteilte Zufallsvariable, dann kann man bei einer Bewegung nach links die Übergangsfunktion von x aus durch U x realisieren, bei einer Bewegung nach rechts durch x + U ( x ). In Ihrem Programm wird also in jedem Schritt eine der beiden Funktionen mit einem Zufallsgenerator mit Wahrscheinlichkeit p = / ausgewählt. a) Welche Rolle spielt die Reihenfolge der Wahl der Zufallsvariablen? Sollte man zunächst alle zufälligen Richtungen simulieren, dann alle U für die Gleichverteilung auf entweder ( , x ) bzw. (x , )? Oder andersherum, erst die Gleichverteilungsvariable U, dann die Richtung? Beides alternierend für jeden Schritt? Die Wahl von U in Abhängigkeit der Richtungswahl? Selbst bei fester Wahl von set.seed werden sich unterschiedliche Folgen von Zahlen ergeben. Die Frage zielt auf die Wahl, die den Prozess am effizientesten (mit der kleinsten Zahl von Aufrufen des Zufallsgenerators) simuliert. b) Diese Markow-Kette hat eine eindeutige stationäre Verteilung mit der Dichte /π x( x) (die Arcussinus-Dichte). Simulieren Sie Realisationen dieser Kette und zeichnen Sie die empirische Dichte der Werte nach Schritten der Markow-Ketten. Vergleichen Sie das Ergebnis mit der Dichtefunktion der stationären Verteilung. Hinweis: density() berechnet einen Kernschätzer, der Befehl curve() kann benutzt werden, die Dichte der stationären Verteilung zu zeichnen. 5) Simulieren Sie die zweidimensionale Markow-Kette, die durch Xt := At Xt + Bt definiert ist. Dabei sollA und BmitWahrscheinlichkeit . dieMatrix Stochastische Modelle bzw. der Vektor A = . . . . , B = + . sein, mit Wahrscheinlichkeit . A = . . . . , B = . . Wählen Sie mindestens Schritte und den Startwert ( . , . ). Zeichnen Sie einen Pfad der Markow-Kette mit dem plot() Befehl und verwenden Sie das Symbol pch=".", also plot(erg,pch="."), wobei erg die Punkte der Kette enthält. 4 10 35 83 _W iS o B eh r P öt te r 2A - B g 4