### Use the mouse to input data:

n = 32
x = rep(0,n)
y = rep(0,n)
plot(c(0,10),c(0,10),type="n",xlab="longitude",ylab="latitude")
for(i in c(1:n)){
loc1 = locator(1)
x[i] = loc1$x
y[i] = loc1$y
points(x[i],y[i])
}





## To save a plot, before making the plot do
##        pdf("myplot.pdf", w = 6, h = 6)
## then after making the plot do
##        dev.off()

plot(c(0,10),c(0,10),type="n",xlab="longitude",ylab="latitude")
points(x,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].

x1 = x
y1 = y

n = length(x) ## n is the number of points
xmax = 10
ymax = 10
a = xmax * ymax ## a = area of observed region.
mu = n/a
alpha = .3
beta =  .1

logl = sum(log(mu+alpha*x+beta*y)) - mu*xmax*ymax -
	alpha*ymax*xmax^2/2 - beta*xmax*ymax^2/2

#### Loglikelihood, for log(lambda) = mu + alpha X + beta Y.
#n = length(x) ## n is the number of points
#xmax = 10
#ymax = 10
#a = xmax * ymax ## a = area of observed region.
#mu = log(n/a)
#alpha = .3
#beta =  .1
#logl = sum(mu+alpha*x+beta*y) - 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.

n = length(x)
mu = n/a

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(mu, .3, .1)
z1 <- nlm(f, pstart, typsize=c(1,.1,.1),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,
#mu = log(n/a)
#
#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(mu, .03, .1)
#z1 <- nlm(f, pstart, typsize=c(1,.1,.1),x=x1,y=y1,print=0)
#pend <- z1$estimate
#pstart <- pend
#z1 <- nlm(f, pstart,x=x1,y=y1,print=0)
#pend <- z1$estimate
##########################################

### PLOT LAMBDA
# 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.

par(mfrow=c(2,2))
image(x2,y2,z1,xlab="x",ylab="y",col=gray((64:20)/64))
mtext(side=3,line=1,"Maximum likelihood estimates")
points(x1,y1)
## Darker means higher rate
######### LEGEND:
zmin <- min(z1)
zmax <- max(z1)
zrng <- zmax - zmin
zmid <- zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid <- seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)

image(x2,y2,z2,xlab="x",ylab="y",col=gray((64:20)/64))
mtext(side=3,line=1,"Original (guessed) parameter estimates")
points(x1,y1)
######### LEGEND:
zmin <- min(z2)
zmax <- max(z2)
zrng <- zmax - zmin
zmid <- zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid <- seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)




### LOAD SM, SPATSTAT, and SPLANCS

## If you use R on the STAT server, the library "splancs" should be
## already installed.  Otherwise you will have to install it.
## try typing 
##      install.packages("sm")
##      library(sm)
##      install.packages("spatstat")
##      library(spatstat)
##      install.packages("splancs")   
##	library(splancs)
## These packages are 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.
## Note: you have to install "sm" before "spatstat"!!! 

library(sm)
library(spatstat)
library(splancs)
## Also if you want you can do
library(help=splancs)      ### briefly describes functions in splancs.
library()		   ### lists available installed libraries.
################################################

b1 <- as.points(x1,y1)
n <- npts(b1)

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

par(mfrow=c(2,2))
bdry <- matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2,byrow=T)
plot(c(0,10),c(0,10),type="n",xlab="x",ylab="y")
lines(bdry)
points(b1,pch="*")
s <- (1:100)*.02
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)