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

### SIMULATE 100 pts distributed uniformly on the square [0,1000] x [0,1000]

xsim <- runif(100)*1000
ysim <- runif(100)*1000

### LOAD IN DATA FROM A TEXT FILE
### (Be sure this file is in the same directory you are in when you start R.)

dat <- scan("datafile.txt")
n <- length(dat)/2
xfile <- dat[1:n]
yfile <- dat[(n+1):(2*n)]

### INPUT DATA DIRECTLY

x1 <- c(752.20, 593.40, 101.40, 357.80, 996.80, 48.05, 355.90, 710.30, 
592.80, 770.10, 442.90, 371.80, 35.16, 901.10, 637.80, 158.60, 765.90, 
326.20, 500.60, 812.20, 441.10, 831.50, 304.40, 743.30, 716.80, 990.70, 
609.70, 825.20, 681.70, 422.40)

y1 <- c(860.90, 447.00, 166.50, 170.90, 730.60, 180.30, 151.10, 277.30, 
514.00, 502.20, 584.80, 370.90, 833.20, 880.20, 447.40, 100.40, 44.70, 
311.90, 837.70, 949.30, 89.24, 405.00, 578.40, 649.50, 27.59, 699.00, 
664.30, 582.40, 590.60, 774.90)

### PLOT DATA

postscript("fig1.ps")
par(mfrow=c(3,2))
plot(xsim,ysim)
plot(xfile,yfile)
plot(x1,y1)
plot(x1,y1,pch="x")
plot(x1,y1,pch="x",cex=.7)
plot(x1,y1,pch="*")
dev.off()

### MORE PLOTTING

postscript("fig2.ps")
par(mfrow=c(2,2))
## A
plot(x1,y1)
points(xsim,ysim,pch="x")
## B
plot(c(0,1000),c(0,1000),type="n",xlab="",ylab="")
points(xsim,ysim,pch="x")
points(x1,y1,pch="o")
mtext(side=1,line=2,"longitude")
mtext(side=2,line=2,"latitude")
mtext(side=3,line=1,cex=2,"Scatterplot")
## C
plot(c(0,1000),c(0,1000),type="n",xlab="",ylab="",axes=F)
lines(c(0,1000,1000,0,0),c(0,0,1000,1000,0))
mtext(side=1,line=-.5,at=0,"119.0",cex=.7)
mtext(side=1,line=-.5,at=1000,"117.0",cex=.7)
mtext(side=2,line=1,at=0,"34.0",cex=.7,las=1)
mtext(side=2,line=1,at=500,"35.0",cex=.7,las=1)
mtext(side=2,line=1,at=1000,"36.0",cex=.7,las=1)
mtext(side=1,line=2,"longitude")
mtext(side=2,line=2,"latitude")
points(xsim,ysim)
## D
plot(c(0,1000),c(0,1000),type="n",xlab="",ylab="",axes=F)
lines(c(-40,1040,1040,-40,-40),c(-40,-40,1040,1040,-40))
points(x1,y1)
points(xsim,ysim,pch="x")
dev.off()

##########################################

### SIMULATE STATIONARY POISSON WITH RATE 1
### on B = [0,10] x [0,50].
par(mfrow=c(1,1))
npts <- rpois(1,(10*50))
xsim2 <- runif(npts)*10
ysim2 <- runif(npts)*50
plot(c(0,10),c(0,50),type="n",xlab="",ylab="")
points(xsim2,ysim2)
mtext(side=1,line=2,"x")
mtext(side=2,line=2,"y")


##########################################

### Compute Loglikelihood, for lambda = mu + alpha X + beta Y
### Given a value of mu, alpha, and beta, and an 
### observation region such as [0,10] x [0,10].
mu <- .00002
alpha <- .00000003
beta <-  .00000001
n <- length(x1)  ## n is the number of points
xmax <- 1000
ymax <- 1000
logl <- sum(log(mu+alpha*x1+beta*y1)) - mu*xmax*ymax -
	alpha*ymax*xmax^2/2 - beta*xmax*ymax^2/2

### Loglikelihood, for log(lambda) = mu + alpha X + beta Y. 
mu <- -10
alpha <- .0001
beta <- .0002
n <- length(x1)  ## n is the number of points
xmax <- 1000
ymax <- 1000
logl <- sum(mu+alpha*x1+beta*y1) - exp(mu)/alpha/beta *
	(exp(alpha*xmax) - 1) * (exp(beta*ymax) - 1)


##########################################

MAXIMUM LIKELIHOOD ESTIMATION:
### first make f(p) the negative loglikelihood,
### where p = (mu, alpha, beta) :
### below is for the case where lambda = mu + alpha X + beta Y.
f <- function(p,x,y){
  if(p[1]+p[2]*xmax+p[3]*ymax < 0) return(99999)
  if(p[1] < 0) return(99999)
  if(p[1]+p[2]*xmax < 0) return(99999)
  if(p[1]+p[3]*ymax < 0) return(99999)
  return(-sum(log(p[1]+p[2]*x+p[3]*y)) + p[1]*xmax*ymax +
	p[2]*ymax*xmax^2/2 + p[3]*xmax*ymax^2/2)
}
pstart <- c(.00002, .00000003, .00000001) 
z1 <- nlm(f, pstart, typsize=pstart,x=x1,y=y1,print=0) 
pend <- z1$estimate 

### TO CHECK, COMPARE THE FOLLOWING:
f(pstart,x1,y1)
f(pend,x1,y1)
### the latter one should be less.
##########################################

### for the case where log(lambda) = mu + alpha X + beta Y,
f <- function(p,x,y){
  -sum(p[1]+p[2]*x1+p[3]*y1) + exp(p[1])/p[2]/p[3] *
	(exp(p[2]*xmax) - 1) * (exp(p[3]*ymax) - 1)
}
pstart <- c(.00002, .00000003, .00000001) 
z1 <- nlm(f, pstart, typsize=pstart,x=x1,y=y1,print=0) 
pend <- z1$estimate 
pstart <- pend
z1 <- nlm(f, pstart, typsize=pstart,x=x1,y=y1,print=0) 
pend <- z1$estimate 
##########################################

### PLOT LAMBDA  
postscript("fig3.ps")
par(mfrow=c(2,2))
plot(x1,y1,xlab="x",ylab="y")
x2 <- c(0:10)/10*xmax
y2 <- c(0:10)/10*ymax
z1 <- matrix(rep(0,(11*11)),ncol=11)
z2 <- z1
for(i in c(1:11)){
for(j in c(1:11)){
z1[i,j] <- pend[1] + pend[2]*x2[i] + pend[3]*y2[j]
z2[i,j] <- pstart[1] + pstart[2]*x2[i] + pstart[3]*y2[j]
}}

##### for log(lambda) = mu + alpha X + beta Y, use exp(z1) 
##### where you see z1 below.  

image(x2,y2,z1,xlab="x",ylab="y",col=gray((0:32)/32))
mtext(side=3,line=1,"Original parameter estimates")
image(x2,y2,z1,xlab="x",ylab="y",col=gray((0:32)/32))
mtext(side=3,line=1,"Original parameter estimates")
points(x1,y1)
image(x2,y2,z2,xlab="x",ylab="y",col=gray((0:32)/32))
mtext(side=3,line=1,"Maximum likelihood estimates")
points(x1,y1)
dev.off()
#########

### PLOT RESIDUALS:  (for lambda = mu + alpha X + beta Y)

postscript("fig4.ps")
yres <- pend[1]*y1 + pend[2]*x1*y1 + pend[3]*y1^2/2
par(mfrow=c(2,2))
plot(x1,y1,xlab="x",ylab="y")
plot(x1,yres,xlab="x",ylab="y")
plot(x1,yres,xlab="x",ylab="y")
yres2 <- pend[1]*ymax + pend[2]*x2*ymax + pend[3]*ymax^2/2
lines(x2,yres2)
b <- max(yres2)
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y")
points(x1,yres)
lines(x2,yres2)
dev.off()

##########################################
### PLOT RESIDUALS:  (for log(lambda) = mu + alpha X + beta Y)
postscript("fig4b.ps")
yres <- exp(pend[1]+pend[2]*x)/pend[3]*(exp(pend[3]*y)-1)
par(mfrow=c(2,2))
plot(x1,y1,xlab="x",ylab="y")
plot(x1,yres,xlab="x",ylab="y")
plot(x1,yres,xlab="x",ylab="y")
yres2 <- exp(pend[1]+pend[2]*x2)/pend[3]*(exp(pend[3]*ymax)-1)
lines(x2,yres2)
b <- max(yres2)
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y")
points(x1,yres)
lines(x2,yres2)
dev.off()


##########################################
### Simulate hom. Poisson with rate 1 in the transformed region.
### (Simulate over rectangle, and thin.)
n2 <- rpois(1,xmax*b)
xsim2 <- runif(n2)*xmax
ysim2 <- runif(n2)*b

keep <- (ysim2 <= pend[1]*ymax + pend[2]*xsim2*ymax + 
	pend[3]*ymax^2/2)

## if you are using log(lambda) = mu + alpha X + beta Y, 
## then replace the above line with: 
## 	keep <- (ysim2 <= exp(pend[1]+pend[2]*xsim2)/pend[3]*
##	(exp(pend[3]*ymax)-1))

postscript("fig5.ps")
par(mfrow=c(2,2))
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y")
points(xsim2,ysim2)
lines(x2,yres2)
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y")
points(xsim2[keep],ysim2[keep])
lines(x2,yres2)
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y")
points(x1,yres)
lines(x2,yres2)
dev.off()


### MISCELLANEOUS

length(x)
y
y[7]
help(plot)
help(par)

### VIEWING AND CONVERTING TO PDF

Quit R, and in Unix do:

ghostview fig1.ps 
ps2pdf fig1.ps