######### After fitting the MLE in rfunctions101.txt,
######### and obtaining parameter estimates "pend".......

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

yres = pend[1]*y1 + pend[2]*x1*y1 + pend[3]*y1^2/2
yres2 = pend[1]*ymax + pend[2]*x2*ymax + pend[3]*ymax^2/2
par(mfrow=c(2,1))
plot(x1,y1,xlab="x",ylab="y",main="original points")
#plot(x1,yres,xlab="x",ylab="y",main="rescaled points")
b <- max(yres2)
plot(c(0,xmax),c(0,b),type="n",xlab="x",ylab="y",main="rescaled points")
points(x1,yres)
lines(x2,yres2)

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






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

library(splancs)

s2 <- (1:100)/100*5  		### 100 numbers between 0 and 5.
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()

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

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

n = length(x1)
d1 = runif(n)
lams = pend[1] + pend[2]*x1 + pend[3]*y1 
bmin = min(pend[1], pend[1] + pend[2]*xmax, pend[1] + pend[3]*ymax, 
pend[1] + pend[2]*xmax + pend[3]*ymax) ## check 4 corners for minimum
keep1 = (d1 < bmin/lams)
thinx1 = x1[keep1]
thiny1 = y1[keep1]
thinn1 = length(thinx1)
par(mfrow=c(2,1))
plot(x1,y1,xlab="x",ylab="y",main="original points")
plot(thinx1,thiny1,xlab="x",ylab="y",main="thinned residuals")


#### THINNED RESIDUALS (for log[lambda] = mu + alpha X + beta Y):

n = length(x1)
d1 = runif(n)
lams = exp(pend[1] + pend[2]*x1 + pend[3]*y1)
bmin1 = min(pend[1], pend[1] + pend[2]*xmax, pend[1] + pend[3]*ymax, 
pend[1] + pend[2]*xmax + pend[3]*ymax) ## check 4 corners for minimum
bmin = exp(bmin)
keep1 = (d1 < bmin/lams)
thinx1 = x1[keep1]
thiny1 = y1[keep1]
thinn1 = length(thinx1)
par(mfrow=c(2,1))
plot(x1,y1,xlab="x",ylab="y",main="original points")
plot(thinx1,thiny1,xlab="x",ylab="y",main="thinned residuals")



### Then can repeat this, and compute the L-function as above....