### CONFIDENCE BOUNDS FOR K-FUNCTION

k4conf <- Kenv.csr(npts(b1), bdry, 1000, s)
plot(c(0,max(s)),c(0,max(k4conf$upper)),
	type="n",xlab="distance",ylab="K4-hat")
points(s,k4,pch="*")
lines(s,k4)
lines(s,k4conf$upper,lty=2)
lines(s,k4conf$lower,lty=2)

L4upper <- sqrt(k4conf$upper/pi) - s
L4lower <- sqrt(k4conf$lower/pi) - s
plot(c(0,max(s)),c(min(L4lower),max(L4upper)),
	type="n",xlab="distance",ylab="L4-hat")
points(s,L4,pch="*")
lines(s,L4)
lines(s,L4upper,lty=2)
lines(s,L4lower,lty=2)

#### Simulate a Neyman-Scott model
library(spatstat)
library(help=spatstat) ## lists all the functions in Spatstat
nclust <-  function(x0, y0, radius, n) {
                                return(runifdisc(n, radius, x0, y0))
                              }
z <- rNeymanScott(10, 0.2, nclust, radius=0.05, n=20)
plot(z,pch=".")

### IMAGE LEGENDS, IN GENERAL:
x1 <- (1:10)
x2 <- (1:5)
x3 <- matrix(runif(50),ncol=5)
## x3 is the data, in this case iid uniform random numbers bet 0 & 1
x4 <- (0:100)/100*(max(x3)-min(x3))+min(x3)
## x4 is a vector of 101 equally-spaced numbers spanning the range of the data.
## For example, if the data values range from exactly 0 to 1,
## then x4 will be {0, .01, .02, .03, ..., .99, 1}
par(mfrow=c(1,2))
image(x1,x2,x3,col=gray((64:0)/64),xlab="longitude",ylab="latitude")
plot(c(0,10),c(.8*min(x4),1.2*max(x4)),type="n",axes=F,xlab="",ylab="")
image(c(-1:1),x4,matrix(rep(x4,2),ncol=101,byrow=T),add=T,col=gray((64:0)/64))
text(2,min(x4),as.character(signif(min(x4),2)),cex=1)
text(2,(max(x4)+min(x4))/2,as.character(signif((max(x4)+min(x4))/2,2)),cex=1)
text(2,max(x4),as.character(signif(max(x4),2)),cex=1)
## I use the function signif() here so that it 
## only prints 2 significantdigits.


###  Kernel smoothing

library(splancs)
b1 <- as.points(x1,y1)
bdry <- matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2,byrow=T)
bdw = 2.2 
z = kernel2d(b1,bdry,bdw)  
## see help(kernel2d)
## b1 should be a "points" dataset, bdry a polygon, bdw = bandwidth


### PLOT KERNEL SMOOTHED DATA WITH THE POINTS, AND WITH A LEGEND
par(mfrow=c(1,2))
image(z,col=gray((64:20)/64),xlab="x",ylab="y")
points(b1)
x4 = (0:100)/100*(max(z$z)-min(z$z))+min(z$z)
plot(c(0,10),c(.8*min(x4),1.2*max(x4)),type="n",axes=F,xlab="",ylab="")
image(c(-1:1),x4,matrix(rep(x4,2),ncol=101,byrow=T),add=T,col=gray((64:20)/64))
text(2,min(x4),as.character(signif(min(x4),2)),cex=1)
text(2,(max(x4)+min(x4))/2,as.character(signif((max(x4)+min(x4))/2,2)),cex=1)
text(2,max(x4),as.character(signif(max(x4),2)),cex=1)
mtext(s=3,l=-3,at=1,"Rate (pts per unit area)")