source("SemiparametricTest_revision_final.R") ##plot MAD of LGCP-B process data0 <- read.csv("./output/output_size_LGCP_bessel_pcf.csv",header=FALSE,sep=",", skipNul = FALSE) for(nu in c(0.5,0.75,1)){ data <- data0[data0[,103]==nu&data0[,104]==0.01&data0[,101]==1,1:100] nrep <- 3 N <- nrow(data) q <- ncol(data) nsim <- N/nrep r <- seq(0,0.15,l=100) gausmod <- RMbessel(var=1, scale=0.01, nu=nu) gfun0 <- gfun_LGCP(r,gausmod) n.r <- length(r) MAD_t <- Mean_t <- array(NA,dim=c(nrep,n.r,2)) #nsim <- 3000 for (i in 1:nrep) { ind <- seq(i,N,by=nrep) data_sub <- (data[ind,])[1:nsim,] dif <- abs(scale(data_sub,center=gfun0,scale=FALSE)) MAD_t[i,,1] <- apply(dif, 2, mean) } data <- data0[data0[,103]==nu&data0[,104]==0.01&data0[,101]==2,1:100] N <- nrow(data) q <- ncol(data) nsim <- N/nrep #nsim <- 3000 for (i in 1:nrep) { ind <- seq(i,N,by=nrep) data_sub <- (data[ind,])[1:nsim,] dif <- abs(scale(data_sub,center=gfun0,scale=FALSE)) MAD_t[i,,2] <- apply(dif, 2, mean) } ylim=c(0,1.1) ylab=expression(MAD[hat(g)](r)) xlab <- "r" if(nu==0.5){ main=expression(bold(paste("(b) LGCP-B (",nu,"=0.5)",sep=''))) filename <- "./fig/Fig3-2-b.pdf" } if(nu==0.75){ main=expression(bold(paste("(c) LGCP-B (",nu,"=0.75)",sep=''))) filename <- "./fig/Fig3-2-c.pdf" } if(nu==1){ main=expression(bold(paste("(d) LGCP-B (",nu,"=1)",sep=''))) filename <- "./fig/Fig3-2-d.pdf" } len1 <- "Constrained PCF estimator (n=1)" len3 <- "Kernel PCF estimator (n=1)" len2 <- "Parametric PCF estimator (n=1)" len4 <- "Constrained PCF estimator (n=2)" len6 <- "Kernel PCF estimator (n=2)" len5 <- "Parametric PCF estimator (n=2)" pdf(filename) par(mai=c(0.9,1,0.5,0.2),mfrow=c(1,1)) plot(r,MAD_t[1,,1],lwd=2,ylim=ylim,xlab=xlab,ylab=ylab,main=main,cex.lab=1.5,cex.main=1.5,cex.axis=1.3,type="l") lines(r,MAD_t[3,,1],lwd=2,lty=1,col="blue") lines(r,MAD_t[1,,2],lwd=2,lty=2,col="black") lines(r,MAD_t[3,,2],lwd=2,lty=2,col="blue") legend(0.05,1.0, legend=c(len1,len3),col=c(1,"blue"),lty=c(1,1),lwd=2,bty = "n") legend(0.05,0.9, legend=c(len4,len6),col=c(1,"blue"),lty=c(2,2),lwd=2,bty = "n") dev.off() } main="(a) True PCFs (LGCP-B)" len1 <- expression(paste("LGCP-B, ",nu,"=0.5")) len2 <- expression(paste("LGCP-B, ",nu,"=0.75")) len3 <- expression(paste("LGCP-B, ",nu,"=1")) ylim <- c(0.8,2.8) ylab=expression(g[0](r)) xlab="r" pdf("./fig/Fig3-2-a.pdf") par(mai=c(0.9,1,0.5,0.2),mfrow=c(1,1)) gausmod <- RMbessel(var=1, scale=0.01, nu=0.5) plot(r, gfun_LGCP(r,gausmod),lwd=2,ylim=ylim,xlab=xlab,ylab=ylab,main=main,cex.lab=1.5,cex.main=1.5,cex.axis=1.3,type="l") gausmod <- RMbessel(var=1, scale=0.01, nu=0.75) lines(r,gfun_LGCP(r,gausmod),type="l",lwd=2,ylim=c(1,3.5),col=2,lty=1) gausmod <- RMbessel(var=1, scale=0.01, nu=1) lines(r,gfun_LGCP(r,gausmod),type="l",lwd=2,ylim=c(1,3.5),col=4,lty=1) abline(h=1) legend("topright", legend=c(len1,len2,len3),col=c(1,2,4),lty=c(1,1,1),lwd=2,bty = "n") dev.off()