#
# Use Shapiro-Wilks W as measure of optimality of two parameter
# Box-Cox transformation
#
boxcox.trans <- function(lam,y) {
    if (lam[1]==0) {
      ynew<-log(y+lam[2])
    }else{
      ynew<-(y+lam[2])^(lam[1])
    } 
    return(ynew)
}
boxcox.grid <- function (y,upper.pow=3,upper.offset=5,step=0.1,plotit=F) {
    if (3>length(y)>5000) {
      print("Error in boxcox: length of y must be between 3 and 5000.")
      return(y)
    }
    require(ctest)
    pow<-c(seq(-upper.pow,-step,by=step),0,seq(step,upper.pow,by=step))
    offset<-c(seq(0,upper.offset,by=step))
    w<-matrix(NA,ncol=length(offset),nrow=length(pow))
    rownames(w)<-pow
    colnames(w)<-offset
    for(i in 1:length(pow)) {
      for(j in 1:length(offset)) {
        w[i,j]<-shapiro.test(boxcox.trans(c(pow[i],offset[j]),y))$statistic
      }
    }
    best.pow<-pow[row(w)[w==max(w,na.rm=T)]][1]
    best.offset<-offset[col(w)[w==max(w,na.rm=T)]][1]
    if (plotit) {
       persp(as.double(rownames(w)),
             as.double(colnames(w)),
             w,ticktype="detailed",
             theta=135,phi=10,col=3,
             xlab="Power",ylab="Offset",zlab="S-W W",shade=0.75)
       title(paste("Peak W at power=",best.pow,", offset=",offset))
    }else{   
       y.new<-boxcox.trans(c(best.pow,best.offset),y)
       return(y.new,lam=c(best.pow,best.offset),w)
    }
}
boxcox.best <- function(y) {
    boxcox.w <- function(lam,y) {
       w<-shapiro.test(boxcox.trans(lam,y))$statistic
       return(w)
    }
    require(ctest)
    res<-optim(par=c(0.0,0.0), fn=boxcox.w, gr=NULL,
               method="Nelder-Mead", lower= -Inf, upper= +Inf,
               control = list(fnscale=-1), hessian = FALSE, y)
    lam<-res$par
    w<-res$value
    ynew<-boxcox.trans(lam,y)
    return(ynew,lam,w)
}    
#
# Looking for optimal transformation for intraclass correlation model
# cluster size 2
#
boxcox2.grid <- function (y1,y2,upper.pow=3,upper.offset=5,step=0.1,plotit=F) {
    if (length(y1)!=length(y2)) {
      print("Error in boxcox2.grid: length of y1 ne y2.")
      return(y1,y2)
    }
    if (3>length(y1)>5000) {
      print("Error in boxcox2.grid: length of y must be between 3 and 5000.")
      return(y1,y2)
    }
    ic1<-c(y1,y2)
    ic2<-c(y2,y1)
    require(ctest)
    pow<-c(seq(-upper.pow,-step,by=step),0,seq(step,upper.pow,by=step))
    offset<-c(seq(0,upper.offset,by=step))
    w<-matrix(NA,ncol=length(offset),nrow=length(pow))
    rownames(w)<-pow
    colnames(w)<-offset
    for(i in 1:length(pow)) {
      for(j in 1:length(offset)) {
        w[i,j]<-shapiro.test(resid(
            lm(boxcox.trans(c(pow[i],offset[j]),ic1) ~
               I(boxcox.trans(c(pow[i],offset[j]),ic2)))
          ))$statistic
      }
    }
    best.pow<-pow[row(w)[w==max(w,na.rm=T)]][1]
    best.offset<-offset[col(w)[w==max(w,na.rm=T)]][1]
    if (plotit) {
       persp(as.double(rownames(w)),
             as.double(colnames(w)),
             w,ticktype="detailed",
             theta=135,phi=10,col=3,
             xlab="Power",ylab="Offset",zlab="S-W W",shade=0.75)
       title(paste("Peak W at power=",best.pow[1],", offset=",best.offset))
    }else{ 
       y1.new<-boxcox.trans(c(best.pow,best.offset),y1)
       y2.new<-boxcox.trans(c(best.pow,best.offset),y2)
       return(y1.new,y2.new,lam=c(best.pow,best.offset),w)
    }
}
boxcox2.best <- function(y1,y2) {
    boxcox2.w <- function(lam,y1,y2) {
       w<-shapiro.test(resid(
            lm(boxcox.trans(lam,y1) ~ I(boxcox.trans(lam,y2)))
          ))$statistic
       return(w)
    }
    ic1<-c(y1,y2)
    ic2<-c(y2,y1)
    require(ctest)
    res<-optim(par=c(0.0,0.0), fn=boxcox2.w, gr=NULL,
               method="Nelder-Mead", lower= -Inf, upper= +Inf,
               control = list(fnscale=-1), hessian = FALSE, ic1, ic2)
    lam<-res$par
    w<-res$value
    y1.new<-boxcox.trans(lam,y1)
    y2.new<-boxcox.trans(lam,y2)
    return(y1.new,y2.new,lam,w)
}