####################################### # # Räumliche Statistik WS 2010 / 2011 # Blatt 3, Aufgabe 1 # ####################################### # Benutze die Prozeduren von Blatt 1####################################### # # Räumliche Statistik WS 2010 / 2011 # Blatt 1, Aufgabe 4 # ##################################### # Erzeuge Realisierungen einer exponentialverteilten # Zufallsvariablen mit Parameter lambda>0, # d.h. Erwartungswert 1 / lambda # Transformation siehe Übungsblatt. getExpRV<-function(lambda) { value<-runif(1) value<- (-1.0 * log(value) / lambda) return(value) } # Erzeuge Realisierungen einer Poissonverteilten # Zufallsvariablen mit Parameter lambda>0, # d.h. Erwartungswert lambda # Algorithmus siehe Übungsblatt. getPoissonRV<-function(lambda) { done<-FALSE value<-0 counter<-0 while(!done) { value<-value + getExpRV(lambda) if(value<=1) { counter<-counter+1 } else { done<-TRUE } } return(counter) } # Nutze die bedingte Gleichverteilungseigenschaft # des Poisson-Prozesses, um eine Realisierung # mit Parameter lambda auf einem rechteckigen # Beobachtungsfenster zu erzeugen. getPoissonPP<-function(window_x_lim, window_y_lim, lambda) { windowSize<-abs(window_x_lim[2]-window_x_lim[1])*abs(window_y_lim[2]-window_y_lim[1]) k<-getPoissonRV(lambda*windowSize) points<-c() for(i in 1:k) { x<-runif(1,min=window_x_lim[1], max=window_x_lim[2]) y<-runif(1,min=window_y_lim[1], max=window_y_lim[2]) points<-rbind(points, c(x,y)) } return(points) } # Visualisiere die Punkte im vorgegebenen Beobachtungsfenster visualizePoints<-function(window_x_lim, window_y_lim, points, filename) { png(filename) plot(window_x_lim, window_y_lim, xlab="", ylab="", main="Realisierung eines Poisson-Prozesses", type="n") par(mfrow=c(1,1), pch=19, lty=3) points(points) dev.off() return(TRUE) } points<-getPoissonPP(c(0,100), c(0,100), 0.01) #print(points) visualizePoints(c(0,100), c(0,100),points, "C:/Lehre/output.png") #source("C:/Lehre/RSBlatt1.r") # Berechne die auf Blatt 2 angebene Teststatistik # für eine Realisierung (points) eines homogenen # Poisson-Prozesses in einem rechteckigen Fenster # und ein vorgegebenes lambda_0>0 calculateTestStatistic<-function(window_x_lim, window_y_lim, points, lambda_0) { windowSize<-abs(window_x_lim[2]-window_x_lim[1])*abs(window_y_lim[2]-window_y_lim[1]) noOfPoints<-dim(points)[1] T<-sqrt(windowSize*windowSize/noOfPoints)*(noOfPoints/windowSize - lambda_0) return(T) } Aufgabe2a<-function(lambda_0) { M<-c(100,200,300,500) # Durchlaufe die verschiedenen Fenstergrößen for(m in M) { statistics<-rep(0,100) # Für jede Fenstergröße finde die Statistiken für # 100 Realisierungen des PPP ... for(i in 1:100) { points<-getPoissonPP(c(0,m), c(0,m), lambda_0) statistics[i]<-calculateTestStatistic(c(0,m), c(0,m), points, lambda_0) } # ... und führe für diese Realisierungen einen Test auf # Standardnormalverteiltheit durch test_result<-shapiro.test(statistics) #print(test_result) # Gib die Testentscheidung aus if(test_result$p.value<0.05) { output<-paste("Die Hypothese der Normalverteiltheit wird beim Fenster mit Seitenlänge",m,"zum Niveau 0.05 verworfen.\n") } else { output<-paste("Die Hypothese der Normalverteiltheit wird beim Fenster mit Seitenlänge ",m," zum Niveau 0.05 nicht verworfen.\n") } cat(output) } } Aufgabe2a(0.001) Aufgabe2b<-function(lambda_0) { Lambda<-c(0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.010) empirical_alpha<-rep(0,length(Lambda)) # Für jedes der vorgegebenen Lambda ... for(l in 1:length(Lambda)) { # ...führe 1000 Tests durch ... for(i in 1:100) { # ... indem zunächst eine Realisierung der Testgröße für lambda bestimmt wird... points<-getPoissonPP(c(0,100), c(0,100), Lambda[l]) t_0<-calculateTestStatistic(c(0,100), c(0,100), points, lambda_0) rank<-1 # ... und dann 99 für lambda_0 for(j in 1:99) { points<-getPoissonPP(c(0,100), c(0,100), lambda_0) t<-abs(calculateTestStatistic(c(0,100), c(0,100), points, lambda_0)) # Bestimme dann den Rang von lambda einfach, indem # gezählt wird, wieviele der 99 Realisierungen für lambda_0 # betragsmäßig kleinere Teststatistiken liefern als lambda if(abs(t)95) { empirical_alpha[l]<-(empirical_alpha[l]+1) } } } # Die empirische Gütefunktion ist also die relative # Häufigkeit des Verwerfens von H_0 empirical_alpha<-empirical_alpha / 100 print(empirical_alpha) png("~/Lehre/RS08/empirical_alpha.png") plot(c(0,0.01), c(0,1), xlab="", ylab="", main="Empirische Gütefunktion", type="n") lines(Lambda,empirical_alpha) dev.off() } Aufgabe2b(0.001) q() Aufgabe2b<-function(lambda_0) { Lambda<-c(0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.010) empirical_alpha<-rep(0,length(Lambda)) # Für jedes der vorgegebenen Lambda ... for(l in 1:length(Lambda)) { # ...führe 1000 Tests durch ... for(i in 1:100) { # ... indem zunächst eine Realisierung der Testgröße für lambda bestimmt wird... points<-getPoissonPP(c(0,1000), c(0,1000), Lambda[l]) t_0<-calculateTestStatistic(c(0,1000), c(0,1000), points, lambda_0) rank<-1 # ... und dann 99 für lambda_0 for(j in 1:99) { points<-getPoissonPP(c(0,1000), c(0,1000), lambda_0) t<-abs(calculateTestStatistic(c(0,1000), c(0,1000), points, lambda_0)) # Bestimme dann den Rang von lambda einfach, indem # gezählt wird, wieviele der 99 Realisierungen für lambda_0 # betragsmäßig kleinere Teststatistiken liefern als lambda if(abs(t)95) { empirical_alpha[l]<-(empirical_alpha[l]+1) } } } # Die empirische Gütefunktion ist also die relative # Häufigkeit des Verwerfens von H_0 empirical_alpha<-empirical_alpha / 100 print(empirical_alpha) png("~/Lehre/RS08/empirical_alpha.png") plot(c(0,0.01), c(0,1), xlab="", ylab="", main="Empirische Gütefunktion", type="n") lines(Lambda,empirical_alpha) dev.off() } Aufgabe2b(0.001) q()