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

xsim <- sort(runif(100)*10)
ysim <- runif(100)*10

### Use the mouse:

n = 22
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])
}


plot(c(0,10),c(0,10),type="n",xlab="longitude",ylab="latitude")
loc1 = locator(n)
points(loc1)
x = loc1$x
y = loc1$y

plot(c(0,10),c(0,10),type="n",xlab="longitude",ylab="latitude")
x[5] = 3.75
points(x,y)


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

sink("datafile.txt")
for(i in 1:n){
cat(x[i])
cat(" ")
cat(y[i])
cat(" ")
}
sink()

write(c(x,y),"datafile2.txt")


dat = scan("datafile2.txt")
n = length(dat)/2
xfile = dat[1:n]
yfile = dat[(n+1):(2*n)]

## alternatively:
dat = scan("datafile.txt")
my.matrix = matrix(dat,ncol=2,byrow=T)
xfile = my.matrix[,1]
yfile = my.matrix[,2]

### INPUT DATA DIRECTLY

x1 = c(7.52, 5.934, 1.014, 3.578, 9.968, 4.805, 3.559, 7.103,
5.928, 7.701, 4.429, 3.718, 3.516, 9.011, 6.378, 1.586, 7.659)

y1 = c(8.609, 4.47, 1.665, 1.709, 7.306, 1.803, 1.511, 2.773,
5.14, 5.022, 5.848, 3.709, 8.332, 8.802, 4.474, 1.004, 4.47)

### PLOT DATA

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


par(mfrow=c(3,2))
plot(x1,y1)
plot(x1,y1,pch="x")
plot(x1,y1,pch="x",cex=.7)
plot(x1,y1,pch="*",cex=1.5)

### MORE PLOTTING

par(mfrow=c(2,2))
## SIMPLEST:
plot(x1,y1)
points(xsim,ysim,pch="x")
## BOTH DATA AND SIMULATION OVERLAID:
plot(c(0,10),c(0,10),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=1.5,"Scatterplot")
## AXES ADJUSTED:
plot(c(0,10),c(0,10),type="n",xlab="",ylab="",axes=F)
lines(c(0,10,10,0,0),c(0,0,10,10,0))
mtext(side=1,line=-.5,at=0,"119.0",cex=.7)
mtext(side=1,line=-.5,at=10,"117.0",cex=.7)
mtext(side=2,line=1,at=0,"34.0",cex=.7,las=1)
mtext(side=2,line=1,at=5,"35.0",cex=.7,las=1)
mtext(side=2,line=1,at=10,"36.0",cex=.7,las=1)
mtext(side=1,line=2,"longitude")
mtext(side=2,line=2,"latitude")
points(xsim,ysim)
## NO AXES AT ALL, and box drawn slightly OUTSIDE range of points
plot(c(0,10),c(0,10),type="n",xlab="",ylab="",axes=F)
lines(c(-.4,10.4,10.4,-.4,-.4),c(-.4,-.4,10.4,10.4,-.4))
points(x1,y1)
points(xsim,ysim,pch="x")

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

### 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].

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 = .00000003
beta =  .00000001

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 = .03
beta =  .01
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, .03, .01)
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"!!! 

## 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

library(splancs)
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)

### 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!!!)

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)



### Other stuff in SPLANCS (getpoly, zoom)

### Note: these don't seem to work all that well.

b5 <- getpoly()  ### use the mouse to draw a polygon.		 
zoom()  	 ### zooms in, but just draws a blank plot with new coordinates.