K = 10 lambda_max <- 0.1 w1 <- 100 w2 <- 100 H <- 11 : 30 lambda_head <- function(X1, X2, x, h, w1, w2){ b1 <- min(x[1] + 0.5 * h, w1) a1 <- max(x[1] - 0.5 * h, 0) b2 <- min(x[2] + 0.5 * h, w2) a2 <- max(x[2] - 0.5 * h, 0) v <- (b1 - a1) * (b2 - a2) cnt <- 0 for(i in 1:length(X1)){ if(X1[i] >= a1 && X1[i] < b1 && X2[i] >= a2 && X2[i] < b2){ cnt <- cnt + 1 } } return(cnt / v) } L <- function(X1, X2, h, w1, w2){ l <- 0 for(i in 1:length(X2)){ B <- (X1 != X1[i]) | (X2 != X2[i]) Y1 <- X1[B] Y2 <- X2[B] l <- l + log(lambda_head(Y1, Y2, c(X1[i], X2[i]), h, w1, w2)) } return(l) } MSE <- function(X1, X2, h, w1, w2){ e <- 0 for(i in 0.5:1:100){ for(j in 0.5:1:100){ e <- e + (lambda_head(X1, X2, c(i,j), h, w1, w2) - i*j*0.00001)^2 } } return(e / (w1 * w2)) } CSV <- matrix(0, K, length(H)) Error <- matrix(0, K, length(H)) for(i in 1:K){ lambda <- lambda_max * w1 * w2 N <- rpois(1, lambda) X1 <- c() X2 <- c() for(j in 1:N){ x1 <- runif(1) * w1 x2 <- runif(1) * w2 if(runif(1)*lambda_max <= x1 * x2 * 0.00001){ X1 <- c(X1, x1) X2 <- c(X2, x2) } } for(k in 1:length(H)){ h = H[k] plot(X1, X2) CSV[i,k] <- L(X1, X2, h, w1, w2) Error[i,k] <- MSE(X1, X2, h, w1, w2) } } print('Estimator for h') for(k in 1: K){ print(10 + which(CSV[k,] == max(CSV[k,]))) } print('Minimal MSE') for(k in 1: K){ print(10 + which(Error[k,] == min(Error[k,]))) }