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

x1 <- sort(runif(100)*10)
y1 <- runif(100)*10

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


n = 28
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
}

x[5] = 3.75
points(x,y)


dat = scan("datafile.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

x11(width = 6, height = 5)
## To get rid of whitespace in between the plots, try:
par(mar = c(5,4,2,2) + .1)
dev.print(dev = postscript, horizontal=F, file="myplot.ps")
## or before making the plot do 
postscript("myplot.ps", width = 6, height = 5)


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

### 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].
n = length(x1) ## 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*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. 
n = length(x1) ## 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*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.
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,alpha,beta) 
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,
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,alpha,beta) 
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  
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((64:20)/64))
mtext(side=3,line=1,"Maximum likelihood estimates")
points(x1,y1)
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)
## Darker means higher rate
#########




### LOAD SPLANCS

## If you use R on the STAT server, the library "splancs" is 
## already installed.  Otherwise you will have to install it.
## try typing install.packages("splancs")  
## 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

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)

### 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(omi=c(1.3,.1,.1,.5))   
par(mfrow=c(3,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)



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