######### 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(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(min(0,k4res,k4resconf$lower),
	max(k4resconf$upper,k4res)),
	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(0,L4reslower,L4res),max(L4resupper,L4res)),
	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)

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

#### 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....
b1 <- as.points(thinx1,thiny1)
k4 <- khat(b1,bdry,s)
L4 <- sqrt(k4/pi)-s
k4conf <- Kenv.csr(thinn1, bdry, 1000, s)
plot(c(0,max(s)),c(0,max(k4conf$upper,k4)),
	type="n",xlab="distance",ylab="K4-hat")
points(s,k4,pch="*")
lines(s,k4)
lines(s,k4conf$upper,lty=2)
lines(s,k4conf$lower,lty=2)

L4upper <- sqrt(k4conf$upper/pi) - s
L4lower <- sqrt(k4conf$lower/pi) - s
plot(c(0,max(s)),c(min(L4lower,L4),max(L4upper,L4)),
	type="n",xlab="distance",ylab="L4-hat")
points(s,L4,pch="*")
lines(s,L4)
lines(s,L4upper,lty=2)
lines(s,L4lower,lty=2)