####################################################################################################
#
#  cdfdr: skew-beta comparison density based nonparametric fdr estimation
#  Date: 2016
#  Cite: Mukhopadhyay, S. (2016) Large-Scale Signal Detection: A Unifying View. Biometrics, 72 325-334.
#
####################################################################################################

library(orthopolynom)
library(fitdistrplus)
library(MASS)

##########################################################################################
#
# The Main Function
#
##########################################################################################

#------------------------------------------------------
CDfdr <- function(z,nb=10){
    #--estimate the empirical null
	qn<-qqnorm(z,plot.it = FALSE)
	emp.par <- as.vector(rlm(qn$y~qn$x,psi = psi.bisquare,method='MM',maxit=100)$coef)
	u <- pnorm(z,mean=emp.par[1],sd=emp.par[2])
	be.coef <- as.numeric(mledist(u,"beta")$estimate)
	uu <- pbeta(u,be.coef[1],be.coef[2])
	S <- Leg.val(uu,nb)
	beta <- apply(S,FUN="mean",2)
	beta <- LP.smooth(beta,length(u),"BIC") 
	d.val <- approxfun(uu,1 + S%*%beta, rule=2,method="constant",f=1 )

    #--p0 estimation----
    dd.val <- function(u){ dbeta(u,be.coef[1],be.coef[2])*d.val(pbeta(u,be.coef[1],be.coef[2]))  }
    D <- dd.val(u)
    nb <- 10
    thres <- seq(1.1,min(3.5,max(D)),length=100)
    T <- rep(NA,length(thres))

    for(k in 1:length(thres)){
		u.thres <- u[which(D<thres[k])]
		S <- Leg.val(u.thres,nb) #
		beta <- LP.smooth(apply(S,FUN="mean",2),length(u),"BIC")
		T[k] <- sum(beta^2)
	}
    s.cut <- thres[min(which(T==min(T)))]
  	signal <- which(D > s.cut)
  	pi0 <- 1 - length(signal)/length(u)
	pi0<- min(1,pi0)
	#--------------------------------------
	
    fdr <- pi0/(dbeta(u,be.coef[1],be.coef[2])*d.val(uu))
    fdr[fdr > 1] <- 1
    F <-list()
    F$fdr <- fdr
    F$fp0 <- c(emp.par,pi0)
    return(F) 

}


##########################################################################################
#
# Two helping functions
#
##########################################################################################

Leg.val <- function(x,d){
     poly <-  slegendre.polynomials(d,normalized=TRUE)
     X <- matrix(NA,length(x),d)
for(j in 1:d){
    X[,j] <- predict(poly[[j+1]],x)
    }
    return(X)
}

LP.smooth <- function(CR,n,method){  
  CR.s <- sort(CR^2,decreasing=TRUE,index=TRUE)$x
  aa <- rep(0,length(CR.s))
   if(method=="AIC"){ penalty <- 2}
   if(method=="BIC"){ penalty <- log(n)}
   aa[1] <- CR.s[1] - penalty/n
   if(aa[1]< 0){ return(rep(0,length(CR)))  }
   if(length(aa)==1){ return(CR) }
   for(i in 2:length(CR.s)){
    aa[i] <- aa[(i-1)] + (CR.s[i] - penalty/n)
    }
  CR[CR^2<CR.s[which(aa==max(aa))]] <- 0
  return(CR)
}