options(object.size=10^100)
###FIRST DEFINE FUNCTION
###doit() fits using local likelihood, window makes the weights, 
###geninv returns the general inverse of a matrix
##THIS IS FUNCTION THAT FITS THE LOCAL-LIKELIHOOD
##The program creates a plot and calls it plot03.ps
doit <- function(){
res <- rep(0,M)
for(t0 in o){
  print(t0)
  if(t0 < (ws+1)/2) { 
    tt <- 1:ws
    ww <- sapply(c((t0-1):0,1:(ws-t0))/(ws-t0+1),function(x) (1-x^3)^3)
    middle <- t0
  } 
    else{ 
      if(t0 > (M - (ws+1)/2 + 1)){
	tt <- (M-ws+1):M
	ww <- sapply(c((t0-(M-ws+1)):1,0:(M-t0))/(t0-(M-ws)),
		     function(x) (1-x^3)^3)
	middle <- ws - M + t0
      }
	else{
	  tt <- (t0 - (ws-1)/2):(t0 + (ws-1)/2)
	  ww <- sapply(c(((ws-1)/2):0,1:((ws-1)/2))/((ws+1)/2),
		       function(x) (1-x^3)^3)
	  middle <- (ws+1)/2
	}
    }
  xx <- rep(Age[tt], N[tt])
  yy <- c()
  for(h in tt)
    yy <- c(yy, rep(1, Y[h]), rep(0, N[h] - Y[h]))
  ww <- rep(ww, N[tt])
  assign("xx",xx,frame=1)
  X <- sapply(0:degree, function(k) xx^k)
  fit1 <- gam(yy ~ X, family = "binomial", control = glm.control(epsilon = 1e-008), weights = ww)
  res[t0] <- fit1$fitted[(cumsum(N[tt]))[middle]]
}
lines(Age,res,lty=LTY)
}


det <- function(M) prod(svd(M, nu = 0, nv = 0)$d)

geninv <- function(M){
  aux <- svd(M)
  auxd <- aux$d
  auxd[auxd!=0] <- 1/auxd
  aux$v%*%diag(auxd,nrow=length(auxd))%*%t(aux$u)
}

Window <- function(n, h=1)
{
        aux <- seq(-0.999, 0.999, len = round(n * h))
        Naux1 <- round(n - length(aux))/2
        Naux2 <- n - length(aux) - Naux1
        c(rep(0, Naux1), (1 - abs(aux)^3)^3, rep(0, Naux2))
}

###READ IN THE DATA ... from loader.dat file
Data <- matrix(scan("loader.dat"),byrow=T,ncol=3)
Age <- Data[,1]
N <- Data[,2]
Y <- Data[,3]
 
###FIND THE CRITERIA FOR EACH WINDOW and EACH model
N0 <- 23
NS <- seq(3,22)
Ks <- c(0:2)
Criteria <- array(0,dim=c(6,length(NS),length(Ks)))
for(i in 1:length(NS)){
  if(NS[i] <= 22) {o <- (N0 - NS[i]):(N0 + NS[i]);M <- length(o)} else {o <- 1:45;M <- 2 * NS[i] + 1}
  xx <- rep(Age[o],N[o])
  yy <- c(); for(h in o) yy <- c(yy,rep(1,Y[h]),rep(0,N[h]-Y[h]))
  ww <- Window(M,1); 
  if(length(ww)>length(o)) ww <- ww[((M-1)/2 - 22):((M-1)/2 + 22)]
  ww <- rep(ww,N[o])
  for(j in 1:length(Ks)){
    X <- sapply(0:Ks[j],function(k) scale01(xx,-1,1)^k)
    fit1 <- glm(yy~X-1,family="binomial",weights=ww)
    p <- fit1$fitted
    AUX <- sweep(X,1,ww,FUN="*")
    AUX <- sweep(AUX,1,p*(1-p),FUN="*")
    AUX0 <- t(XX)%*%XX
    AUX1 <- t(AUX)%*%X
    AUX <- sweep(AUX,1,ww,FUN="*")
    AUX2 <- t(AUX)%*%X 
    NP <- sum(diag(geninv(AUX1) %*% AUX2))	
    l0 <- sum(ww*(yy*log(p) + (1-yy)*log(1-p)))
    Criteria[1,i,j] <- (-2*l0 + 2*j)/length(N)
    Criteria[2,i,j] <- (-2*l0 + 2*NP)/sum(ww)
    Criteria[3,i,j] <- (-2*l0 + j*log(det(AUX0)))/length(N)
    Criteria[4,i,j] <- (-2*l0 + log(det(AUX1)))/sum(ww)
    Criteria[5,i,j] <- (-2*l0 + j*log(det(AUX0)) + 2*j)/length(N)
    Criteria[6,i,j] <- (-2*l0 + log(det(AUX1)) + 2*NP)/sum(ww)
  }
}

###MAKE THE PLOTS of CRITERIA VALUES
postscript("plot03.ps",horizontal=T)
Title <- c("WAIC","WBIC","WCAICF")
par(mfrow=c(2,2), mar=c(3.5,2,2,.2),mgp=c(2,.5,0))
for(i in 1:3){
  plot(2*NS+1,Criteria[i*2,,1],ylim=range(Criteria[i*2,,]),type="l",main=Title[i],yaxt="n",xlab="Window Size in Years")
  lines(2*NS+1,Criteria[i*2,,2],lty=2)
  lines(2*NS+1,Criteria[i*2,,3],lty=3)
  if(i==1) legend(6,.713,legend=c("constant","line","parabola"),lty=c(1,2,3),cex=.8)
}

###FIND WHAT WINDOW SIZE, MODEL MINIMIZE THE CRITERIA
aux <- matrix(0,6,2)
for(k in 1:6){
  x <- order(as.vector(Criteria[k,,]))[1]
  i <- ceiling(x/length(NS))
  aux[k,]<-c(i,NS[x - (i-1)*length(NS)])
}

###MAKE PLOT OF FITTED VALUES
M <- length(N)
o <- 1:M
plot(Age,Y/N,xlab="Age in Years",ylab="Percent Survived")
for(i in c(4,6)){
  degree <- aux[i,1]-1;ws <- 2*aux[i,2]+1;LTY <- (i-2)/2
  doit()
}
degree <- 0; ws <- 13;LTY <- 3##these are local minima
doit()
degree <- 1; ws <- 15;LTY <- 4
doit()
legend(55,.92,legend=c("WAIC,WBIC","WCAICF","constant with ws=13","line with ws=15"),lty=c(1,2,3,4),cex=.8)	
 		      
dev.off()