####################################### # # Raeumliche Statistik WS 2010 / 2011 # Blatt 8, Aufgabe 1 # ####################################### # Benutze die Prozeduren von Blatt 1 & Blatt 3 & Blatt 4 & Blatt 7 # Funktion zur Berechnung von Abständen distance<-function(x_1, y_1, x_2, y_2) { d<-sqrt((x_1-x_2)^2 + (y_1 - y_2)^2) return(d) } # 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) } # Erzeuge eine Realisierung eines Matern-Cluster-Prozesses im vorgegebenen # Fenster mit Elternintensitaet lambda_0, Clusterintensitaet lambda_1, und # Clusterradius R getMaternClusterPP<-function(window_x_lim, window_y_lim, lambda_0, lambda_1, R) { # vergroessere das Fenster in jede Richtung um Radius R, um # keine Randeffekte entstehen zu lassen. plusWindow_x_lim<-c(window_x_lim[1]-R, window_x_lim[2]+R) plusWindow_y_lim<-c(window_y_lim[1]-R, window_y_lim[2]+R) # Simuliere den Eltern-Prozess parentPoints<-getPoissonPP(plusWindow_x_lim,plusWindow_y_lim, lambda_0) allPoints<-c() # Fuer jeden Elternpunkt... for(i in 1:(dim(parentPoints)[1])) { # ... bestimme das ein Fenster um den Elternpunkt mit Kantenlänge 2R ... localWindow_x_lim<-c(parentPoints[i,1]-R , parentPoints[i,1]+R) localWindow_y_lim<-c(parentPoints[i,2]-R , parentPoints[i,2]+R) # ... und simuliere in diesem Fenster den Cluster. daughterPoints<-getPoissonPP(localWindow_x_lim, localWindow_y_lim, lambda_1) # Die Tochterpunkte werden nur dann als Teil der Realisierung des Cluster-Prozesses # uebernommen, wenn sie sowohl naeher als R am Elternprozess liegen als auch innerhalb # des urspruenglich vorgegebenen Fensters. for(j in 1:(dim(daughterPoints)[1])) { if(distance(daughterPoints[j,1],daughterPoints[j,2],parentPoints[i,1],parentPoints[i,2])=window_x_lim[1]) && (daughterPoints[j,1]<=window_x_lim[2]) && (daughterPoints[j,2]>=window_y_lim[1]) && (daughterPoints[j,2]<=window_y_lim[2])) { allPoints<-rbind(allPoints, daughterPoints[j,]) } } } } return(allPoints) } ######################## # a) ######################## # Öffne 3 Grafik-Fenster in einer Zeile par(mfrow=c(1,3)) # simuliere die angegebenen Matern-Cluster-Prozesse lambda0 <- c(1E-5, 1E-4, 1E-3) for(i in 1:3){ lambda1 <- 0.001 / (50^2*pi*lambda0[i]) clusterPoints<-getMaternClusterPP(c(0,1000),c(0,1000), lambda0[i], lambda1, 50) plot(clusterPoints) } ####################### # b) ####################### # Prozedur zur randkorrigierten Schätzung der # nächsten Nachbar-Abstands-Verteilungsfunktion getNearestNeighborDistribution <- function(pointProcess, window_x_lim, window_y_lim, bandWidth, startValue){ # berechne die nächsten Nachbarn # sowie Abstände zum Rand nearestNeighborDistance <- c() boundaryDistance <- c() for (i in 1:length(pointProcess[,1])){ nearestNeighborDistance[i] <- getNearestNeighbor(pointProcess, i) boundaryDistance[i] <- getBoundaryDistance(window_x_lim, window_y_lim, pointProcess[i,1], pointProcess[i,2]) } # berechne zunächst die Normierungskonstante lambda_H <- 0 for (i in 1:length(pointProcess[,1])){ if(nearestNeighborDistance[i] < boundaryDistance[i]){ lambda_H <- lambda_H + 1/((window_x_lim[2]-window_x_lim[1]-2*nearestNeighborDistance[i])*(window_y_lim[2]-window_y_lim[1]-2*nearestNeighborDistance[i])) } } #berechne jetzt für jeden Radius die NNAVF nearestNeighborDistribution <- c() localNearestNeighbor <- 0 radius <- startValue while (localNearestNeighbor < 1-0.01){ localNearestNeighbor <- 0 for (i in 1:length(pointProcess[,1])){ if(nearestNeighborDistance[i] < boundaryDistance[i] & nearestNeighborDistance[i] <= radius){ localNearestNeighbor <- localNearestNeighbor + 1 / ((window_x_lim[2]-window_x_lim[1]-2*nearestNeighborDistance[i])*(window_y_lim[2]-window_y_lim[1]-2*nearestNeighborDistance[i])) } } localNearestNeighbor <- localNearestNeighbor / lambda_H nearestNeighborDistribution <- rbind(nearestNeighborDistribution, c(radius, localNearestNeighbor)) radius <- radius + bandWidth } # for (i in 1:length(nearestNeighborDistance)){ # print(nearestNeighborDistance[i]) # } return(nearestNeighborDistribution) } # Berecnet den nächsten Nachbarn des i-ten Punktes # zu den anderen Punkten in pointProcess, getNearestNeighbor <- function(pointProcess, i){ # Anfangswert auf unendlich setzen smallestDistance <- Inf for (k in 1:length(pointProcess[,1])){ if (k != i){ localDistance <- distance(pointProcess[k,1], pointProcess[k,2], pointProcess[i,1], pointProcess[i,2]) # print(pointProcess[k,1]) # print(pointProcess[k,2]) # print(pointProcess[i,1]) # print(pointProcess[i,2]) # print(localDistance) if ( localDistance < smallestDistance){ smallestDistance <- localDistance } } } return(smallestDistance) } getBoundaryDistance <- function(window_x_lim, window_y_lim, x,y){ # berechne für alle 4 Seiten des Fensters den Abstand distances <- 0 distances[1] <- distance(window_x_lim[1],y,x,y) distances[2] <- distance(window_x_lim[2],y,x,y) distances[3] <- distance(x,window_y_lim[1],x,y) distances[4] <- distance(x,window_y_lim[2],x,y) return(min(distances)) } #################### # c) #################### # Öffne 3 Grafik-Fenster in einer Zeile par(mfrow=c(1,3)) # simuliere die angegebenen Matern-Cluster-Prozesse lambda0 <- c(1E-5, 1E-4, 1E-3) for(i in 1:3){ lambda1 <- 0.001 / (50^2*pi*lambda0[i]) clusterPoints<-getMaternClusterPP(c(0,1000),c(0,1000), lambda0[i], lambda1, 50) plot(function(x) 1-exp(-0.001*pi*x*x), xlim=c(0,50),col="red") points(getNearestNeighborDistribution(clusterPoints, c(0,1000), c(0,1000), 1, 0),type="l",add=T) }