The following commands worked on 5/30/01, in R for Unix version 1.2.1 
which is the version currently installed on the STAT server. 

################################################
### PLOT SIZES

par(omi=c(1.3,.1,.1,.5))   ### Fits the plots in the region so it's
			   ### not cut off on the right.
################################################
### LOAD SPLANCS

## If you use R on the STAT server, the library "splancs" is 
## already installed.  Otherwise you will have to install it.
## splancs is available at http://cran.r-project.org/ , 
## click on "contributed packages" which is at the top of the page, 
## under "ALL PLATFORMS", then look alphabetically for "splancs".

## Within R, to load the splancs library, type
library(splancs)
## Also if you want you can do
library(help=splancs)      ### briefly describes functions in splancs.
library()		   ### lists available installed libraries.
################################################

################################################
### MISC STUFF IN SPLANCS

## as.points, csr, Ghat, npts, pdense

postscript("fig6.ps")
par(omi=c(1.3,.1,.1,.5))   
par(mfrow=c(2,1))
b1 <- as.points(x1,y1)

### Simulate CSR in a polygon:
b2 <- matrix(c(0,0,1000,0,1000,1000,0,800,500,500,0,0),ncol=2,byrow=T)
b3 <- csr(b2,50)
plot(c(0,1000),c(0,1000),type="n",xlab="x",ylab="y")
points(b3,pch="*")
lines(b2)

### G-hat (s) = proportion of pts within s of another point.
s <- (1:40)*5
b4 <- Ghat(b3,s)
plot(s,b4,pch="*")
lines(s,b4)
dev.off()

### Count number of points in the dataset
n <- npts(b1)

### Calculate the avg number of pts per unit area in a polygon.
pdense(b3,b2)	 

################################################
### K FUNCTIONS
### (Remember to LOAD SPLANCS first!!!)

postscript("fig7.ps")
par(omi=c(1.3,.1,.1,.5))   
par(mfrow=c(3,2))
bdry <- matrix(c(0,0,1000,0,1000,1000,0,1000,0,0),ncol=2,byrow=T)
plot(c(0,1000),c(0,1000),type="n",xlab="x",ylab="y")
lines(bdry)
points(b1,pch="*")
s <- (1:100)*2  ### this is the same as:    s <- seq(2,200,2)
k4 <- khat(b1,bdry,s)
plot(s,k4,xlab="distance",ylab="K4-hat",pch="*")
lines(s,k4) 
L4 <- sqrt(k4/pi)-s
plot(s,L4,xlab="distance",ylab="L4-hat",pch="*")
lines(s,L4)

### 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)
dev.off()

################################################
### K FUNCTIONS & CONF BOUNDS ON RESIDUALS

### 1st, find a polygon similar to boundary of transformed region.
### we previously had (x2, yres2) defining the boundary.

bdry2 <- cbind(c(0,x2,xmax,0),c(0,yres2,0,0))

### show the boundary to make sure it's right
postscript("fig8.ps")
par(omi=c(1.3,.1,.1,.5))   
par(mfrow=c(3,2))
plot(c(0,xmax),c(0,max(yres2)),type="n",xlab="x",ylab="y")
points(bdry2,pch="*")
lines(bdry2)
points(x1,yres,pch="x")  	### shows the resids

s2 <- (1:100)/100*50  		### 100 numbers between 0 and 50.
resids <- as.points(x1,yres)
k4res <- khat(resids,bdry2,s2)
plot(s2,k4res,xlab="distance",ylab="K4-hat of residuals",pch="*")
lines(s2,k4res) 
L4res <- sqrt(k4res/pi)-s2
plot(s2,L4res,xlab="distance",ylab="L4-hat of residuals",pch="*")
lines(s2,L4res)

### CONFIDENCE BOUNDS FOR K-FUNCTION

k4resconf <- Kenv.csr(npts(resids), bdry2, 1000, s2)
plot(c(0,max(s2)),c(0,max(k4resconf$upper)),
	type="n",xlab="distance",ylab="K4-hat of residuals")
points(s2,k4res,pch="*")
lines(s2,k4res)
lines(s2,k4resconf$upper,lty=2)
lines(s2,k4resconf$lower,lty=2)

L4resupper <- sqrt(k4resconf$upper/pi) - s2
L4reslower <- sqrt(k4resconf$lower/pi) - s2
plot(c(0,max(s2)),c(min(L4reslower),max(L4resupper)),
	type="n",xlab="distance",ylab="L4-hat of residuals")
points(s2,L4res,pch="*")
lines(s2,L4res)
lines(s2,L4resupper,lty=2)
lines(s2,L4reslower,lty=2)
dev.off()

################################################
### LOGLIKELIHOOD, AIC, ETC.

loglike <- -1*f(pend,x1,y1)  ### this is loglikelihood 
aic <- -2*loglike + 2*length(pend) ### this is the AIC.
## You want the loglikelihood as large as possible, 
## and the AIC as small as possible.

## Note that, for the case where lambda = mu + alpha X + beta Y,
## the following terms in the calculation of f: 
## mu*xmax*ymax + alpha*ymax*xmax^2/2 + beta*xmax*ymax^2/2
## should SUM to approximately  n  for the correct choice of parameters.


## ( Remember:  p[1] = mu, p[2] = alpha, p[3] = beta. )

## This fact can help you get initial starting guesses for 
## mu, alpha, and beta.
## Further, the term mu*xmax*ymax is the expected number of observed
## points just due to the background rate, and the next two terms 
## are the expected additional numbers of points due to trends in the
## x and y directions, respectively.  You can thus look at the 
## plot of your points to get initial guesses for mu, alpha, and beta,
## and calculate mu*xmax*ymax, alpha*ymax*xmax^2/2, & beta*xmax*ymax^2/2
## to make sure these guesses seem about right.

## For the model log(lambda) = mu + + alpha X + beta Y,
## the PRODUCT   
## exp(mu)/alpha/beta * (exp(alpha*xmax) - 1) * (exp(beta*ymax) - 1)
## should be about   n   , if mu, alpha, and beta are guessed right.
## Taking logs of both sides, you get:
## 	[ mu ]  
## 	+ [ log{exp(alpha*xmax)-1} - log(alpha) ] 
## 	+ [ log{exp(beta*ymax) - 1} - log(beta) ]
## should add up to about log(n).
## This first term, [mu], is the contribution of the background rate
## to log(n), the second term is the contribution of the trend in x,
## and the third term is the contribution of the trend in y.

################################################
### Cool but irrelevant plot stuff in SPLANCS (getpoly, zoom)

b5 <- getpoly()  ### you use the mouse to draw a polygon.  
		 ### Cool, but doesn't seem to work that well. 

zoom()  	 ### zoom in.