### G1 # Schreiben Sie eine Funktion mit Namen "fib", die bei Eingabe einer natürlichen Zahl n # die ersten n Fibonaccizahlen # berechnet, speichert und plottet, wobei die Fibonaccizahlen eine rekursive Folge mit # a_0=0, a_1=1 und a_n=a_{n-1}+a_{n-2} sind # Hinweis: # Nutzen Sie dazu ein Leerobjekt, das Sie befüllen. ### Standardbeispiel für Monte-Carlo: # Berechnung der Kreiszahl pi: # Erinnerung: r^2 pi ist die Fl?che eines Kreises, also 1/4 pi die Fl?che des Viertelkreises in einem Einheitsquadrat # Erzeugt man nun eine Anzahl n von gleichverteilten Punkten in diesem Quadrat, so kann man mit # Vergr??erung dieser Anzahl pi immer n?her als den vierfachen Anteilswert der Punkte im Viertelkreis bestimmen. # Vorgehensweise: # Zufallsgenerator setzen: set.seed(1234) # Anzahl der gleichverteilten Punkte n <- 100000 # Ben?tigen Punktepaare (x,y) im Einheitsquadrat mit L?nge = Breite = 1 x <- runif(n, min=0, max=1) y <- runif(n, min=0, max=1) plot(x,y) # Einzeichnen des Viertelkreises: # Aus x^2+y^2=1 folgt f(x)=sqrt(1-x^2) z <- seq(0,1,0.01) lines(z, sqrt(1-z^2), col="red") # Berechnung des Anteils der Punktepaare, bei denen x^2+y^2 <=1 counter <- x^2+y^2<=1 inCircle <- sum(counter) portion <- inCircle / n pi <- 4* portion pi # wahres pi: 3.14159265359... # Verbesserung: schreibe Funktion, die in Abh?ngigkeit der verteilten Punktezahl n pi wiedergibt: calcAndGraphPi <- function(n) { x <- runif(n, min=0, max=1) y <- runif(n, min=0, max=1) plot(x,y) # Einzeichnen des Viertelkreises: # Aus x^2+y^2=1 folgt f(x)=sqrt(1-x^2) z <- seq(0,1,0.01) lines(z, sqrt(1-z^2), col="red") # Berechnung des Anteils der Punktepaare, bei denen x^2+y^2 <=1 portion <- table(x^2+y^2<=1)["TRUE"] pi <- 4*portion/n return(pi) } piEst <- calcAndGraphPi(2000) # Auch wichtig: Wie h?ngt der Wert pi von Anzahl der verteilten Punkte ab? calcPi <- function(n) { x <- runif(n, min=0, max=1) y <- runif(n, min=0, max=1) # Berechnung des Anteils der Punktepaare, bei denen x^2+y^2 <=1 portion <- table(x^2+y^2<=1)["TRUE"] pi <- 4*portion/n return(pi) } n <- matrix(seq(10000,100000, 10000), ncol=1) points <- apply(n, MARGIN=1, FUN=calcPi) plot(n, points) abline(h=3.14159265359) mean(points) -3.14159265359 # vgl: Schwaches Gesetz der gro?en Zahlen ################################################# ### Simulation zur Evaluation von Sch?tzeigenschaften ################################################# # # Beispiel Leeb und P?tscher (2005) "Model Selection: Facts and Fiction" http://journals.cambridge.org/action/displayFulltext?type=1&fid=280990&jid=ECT&volumeId=21&issueId=01&aid=280988 # Frage: Welchen Einfluss hat Modellwahl auf die Sch?tzeigenschaften # im Regressionsmodell mit 2 Regressoren? # # DGP: y = alpha * x1 + beta * x2 + u # Man sei nur an alpha interessiert und muss entscheiden, ob x2 ins Modell geh?rt # oder nicht! # Zielkonflikt: # x2 Weglassen -> kleine Sch?tzvarianz für alpha, aber evtl. omitted variable bias! # Idee: Sch?tzen von beiden Modellen für VIELE Beobachtungen mit FESTEN X, VARIABLEN U # Simulieren Sie für folgende (feste!) X alpha <- 1 beta <- .5 reps <- 1000 # Anzahl Replikationen n <- 80 # Stichprobengr??e rho <- 0.8 # Korrelation zwischen x1 und x2 sig <- c(sdx1=1,sdx2=1,sdu=1) # Stds von x1, x2 und u # Berechne Kovarianzmatrix von (x1, x2) = 2x2 Matrix mit den vier Eintr?gen (CovX <- matrix( c( sig[1]^2, sig[1]*sig[2]*rho, sig[1]*sig[2]*rho, sig[2]^2), ncol = 2)) # EINE (!) n x 2 Matrix X wird ereugt, die sofort mit der Wurzel der Korrelationsstruktur, # also der Choleski Zerlegung der Varianz-Kovarianzmatrix multipliziert wird (X <- matrix(rnorm(2*n),n)%*%t(chol(CovX))) # Initialisierung der Vektoren für die gesch?tzten alphas alpha_u <- numeric(reps) # unrestringiertes Modell mit zwei Regressoren alpha_r <- numeric(reps) # restringiertes Modell mit einem Regressor alpha_sc <- numeric(reps) # Vorher mit Schwarz Kriterium gew?hltes Modell # eigentliche Simulation: # der Prozess wird nun reps mal wiederholt, wobei die variablen u und somit auch variablen y # jedes mal neu zu Beginn eines Schleifendurchlauf erzeugt werden set.seed(1) for (i in 1:reps){ # Simulation von u, y und Erstellen eines Datensatzes (alpha=1, beta=1 !) u <- rnorm(n,sd=sig[3]) # Datensatz hat den Vorteil, das man ihn bei verschiedenen Modellen sch?ner ansprechen kann data <- data.frame(y=alpha*X[,1]+beta*X[,2]+u,x1=X[,1],x2=X[,2]) # Sch?tzen des unrestringierten und restringierten Modells U <- lm(y~-1+x1+x2,data=data) R <- lm(y~-1+x1,data=data) # Koeffizientensch?tzer für alpha im jeweiligen Modell alpha_u[i] <- U$coef[1] alpha_r[i] <- R$coef[1] # Modellwahl mit SC und Speichern des resultierenden Sch?tzers # Falls der BIC Wert des unrestringierten Modells kleiner ist, als der des restringierten Modells, # soll die unrestringierte sch?tzung in alpha_sc gespeichert werden, ansonsten die andere if (AIC(U,k=log(n)) "Verteilungstest") # Frage: War die Notenverteilung der Methoden-Lernzielkontrolle gleichverteilt? n <- 37 # Teilnehmer k <- 13 # Notenklassen x <- c(1,7,3,3,2,3,3,5,2,4,2,2,0) # Beobachtete Notenverteilung x0 <- rep(37/13,13) # Theoretische Notenverteilung bei Gleichverteilung x0 (test_stat <- sum((x-x0)^2/x0)) # Teststatistik des Xi^2 Anpassungstests (pchisq(test_stat,df=k,lower.tail=FALSE)) # P-Value (asymptotisch) # Aber: Xi^2 Verteilung gilt nur asymptotisch, die Verteilung unter # H_0 l?sst sich aber exakt bestimmen: Per Simulation! # Aufgabe nun: Simulieren von 9999 Realisationen der Teststatistik bei Gültigkeit von H0. # Bestimmen des kritischen Werts für einen Test von H0: ### Xi^2 Verteilungstest: Verwenden der Exakten Stichprobenverteilung unter H0: # Anzahl Replikationen Reps <- 9999 # Initialisieren: Vektor der Simulierten Teststatistiken test_stat_sim <- numeric(Reps) for (i in 1:Reps) { # Ziehe zuf?llig Noten aus Gleichverteilung noten <- c(1,1.3,1.7,2.0,2.3,2.7,3.0,3.3,3.7,4.0,4.3,4.7,5.0) noten.stud <- sample(noten,37,replace=TRUE) # Berechne H?ufigkeiten x_sim <- rep(NA,13) for (j in 1:13) x_sim[j] <- sum(noten.stud==noten[j]) # Simulierte Teststatistik test_stat_sim[i] <- sum((x_sim-x0)^2/x0) } # Simulierter P-Value: Anteil der Simulierten Teststatistiken, die gr??er sind # als die beobachtete mean(test_stat_sim>test_stat) # Ist die Chi^2-Verteilung eine gute Approximation? plot(density(test_stat_sim), main="Simulierte Dichte der Teststatistik unter H0 und Xi^2 Approximation") lines(seq(0,max(test_stat_sim),0.1),dchisq(seq(0,max(test_stat_sim),0.1),df=k),col="orange") ### # Weitere M?glichkeiten: Zeitreihenprozesse wie in ?ko2 zur graphischen veranschaulichung ### Aufgabe G2 # Verwenden Sie den datengenerierenden Prozess zur Modellselektion von vorher. # Führen Sie anstelle des Modellvergleichs mit SC in jeder Iteration # einen statistischen Test auf beta=0 durch und w?hlen Sie das restringierte # Modell, wenn dieser Test abgelehnt wird. Vergleichen Sie das Ergebnis mit dem von vorher. # Zeichnen Sie auch diese Verteilung wieder ein. ### Aufgabe G3 # Führen Sie nun jeweils mit dem gew?hlten Modell einen Test auf H_0: alpha=1 # durch. Führen Sie den gleichen Test im unrestringierten Modell durch. # Veranschaulichen Sie das Ergebnis.