Z<-50
r<-seq(0,0.25,0.01)
for (z in 1:Z) {
par(mfrow=c(2,2))
kfun<-rep(0,length(r))
err<-rep(0, length(r))
n<-50
x<-runif(n)
y<-runif(n)
xx<-c(x,x-1,x-1,x-1,x,x,x+1,x+1,x+1)
yy<-c(y,y-1,y,y+1,y-1,y+1,y-1,y,y+1)
plot(xx,yy)
xxx<-rep(xx,n*9)
yyy<-rep(yy,n*9)
dim(xxx)<-c(n*9,n*9)
dim(yyy)<-c(n*9,n*9)
txxx<-t(xxx)
tyyy<-t(yyy)
d<-sqrt((xxx-txxx)^2+(yyy-tyyy)^2)
dd<-d[1:n, ]
for (i in 1:length(r)){
rr<-dd<=r[i]
kfun[i]<-sum(rr)/(n^2)


library(spatstat)
Pt<-ppp(x,y)
K<-Kest(Pt,r)
K1<-K$border
K2<-K$iso
K3<-K$Ripley
K4<-K$trans


newd<-pairdist(Pt)
diag(newd)<- Inf
X<-rep(x,n)
Y<-rep(y,n)
dim(X)<-c(n,n)
dim(Y)<-c(n,n)

ones<-matrix(rep(1,n^2), nrow=n)
aa1<- t(X)>=newd & t(Y)>=newd & t(X)<=(ones-newd) & t(Y)<=(ones-newd)
aa2<- newd<=r[i]
vol<-(ones-2*newd)*(ones-2*newd)
V<-1/vol
kk<-V %*% aa1 %*% aa2
newk[i]<-sum(kk)/(n^2)
}

err1[z, ]<-(K1-pi*r^2)^2
err2[z, ]<-(K2-pi*r^2)^2
err3[z, ]<-(K3-pi*r^2)^2
err4[z, ]<-(K4-pi*r^2)^2
err5[z, ]<-(kfun-pi*r^2)^2
err6[z, ]<-(newk-pi*r^2)^2

rmse1<-sqrt((sum(err1[z, ]))/z)
rmse2<-sqrt((sum(err2[z, ]))/z)
rmse3<-sqrt((sum(err3[z, ]))/z)
rmse4<-sqrt((sum(err4[z, ]))/z)
rmse5<-sqrt((sum(err5[z, ]))/z)
rmse6<-sqrt((sum(err6[z, ]))/z)

}