## ALGORITHM FOR FINDING THE MLE OF A U-SHAPED HAZARD FUNCTION
## August 30, 2008
## Modified November 11, 2008

## THIS FUNCTION FINDS THE CONSTRAINED MLE

find.c.mle	<-	function(x, min.int, plot=FALSE){

	x	<-	sort(x)
	n	<-	length(x)
	F	<- 	ecdf(x)
	S	<-	function(t) 1-F(t)
	AA	<-	cumsum(n*sapply(c(0, x[1:(n-1)]), S)*c(x[1],diff(x)))
	AA	<-	sort(unique(AA)) #NEW
	n	<- 	length(AA) #NEW
	min.int	<-	min.int
	
	if((min.int==1)&&(plot==TRUE)){
		plot(cbind(0,0), cex = 1, pch=21, bg="black", ylab="", xlab="", xlim=c(0,n), ylim=c(0, max(AA)))
		points(cbind(0,0), pch=21, cex=2)
		}

	# find the decreasing bit
	if(min.int>1){
		dec.points	<-	cbind(c(0, 1:(min.int-1),0),c(0, AA[1:(min.int-1)], AA[min.int-1]))

		hpts <- chull(dec.points)
		hpts <- c(hpts, hpts[1])
		if(plot==TRUE){
			plot(dec.points, cex = 1, pch=21, bg="black", ylab="", xlab="", xlim=c(0,n), ylim=c(0, AA[n]))
			lines(dec.points[hpts, ])
			points(dec.points[hpts, ], pch=21, cex=2)
			}

		hpairs	<-	matrix(0,length(hpts)-1,2)
		for(i in 1:(length(hpts)-1)){
			hpairs[i,]	<-	c(hpts[i+1],hpts[i])
			}
		m		<-	length(dec.points[,1])
		hpairs	<-	hpairs[which(hpairs[,1]!=m),]
		hpairs	<-	hpairs[which(hpairs[,2]!=m),]
		if(length(hpairs)==2){
			hpairs<-	matrix(hpairs,1,2)
			} else 	{
					hpairs	<-	hpairs[order(hpairs[,1]),]
					}
				
		dec.slopes	<-	(hpairs[,2]-hpairs[,1])/(dec.points[hpairs[,2],2]-dec.points[hpairs[,1],2])

		dec.mle	<-	numeric()

		if(length(dec.slopes)>0){
			for (i in 1:length(dec.slopes)){
				dec.mle	<-	c(dec.mle, rep(dec.slopes[i], dec.points[hpairs[i,2],1]-dec.points[hpairs[i,1],1]))
				}
			}

		} else {dec.mle <- numeric()}
		

	if((min.int==n)&&(plot==TRUE)){
		points(cbind(n,AA[n]), cex = 1, pch=21, bg="black")
		points(cbind(n,AA[n]), pch=21, cex=2)
		}


	# find the increasing bit
	if(min.int<n){
		inc.points	<-	cbind(c(min.int:n, n),c(AA[min.int:n],AA[min.int]))

		hpts <- chull(inc.points)
		hpts <- c(hpts, hpts[1])
		if(plot==TRUE){
			points(inc.points, cex = 1, pch=21, bg="black")
			lines(inc.points[hpts, ])
			points(inc.points[hpts, ], pch=21, cex=2)
			}

		hpairs	<-	matrix(0,length(hpts)-1,2)
		for(i in 1:(length(hpts)-1)){
			hpairs[i,]	<-	c(hpts[i],hpts[i+1])
			}
		m		<-	length(inc.points[,1])
		hpairs	<-	hpairs[which(hpairs[,1]!=m),]
		hpairs	<-	hpairs[which(hpairs[,2]!=m),]
		if(length(hpairs)==2){
			hpairs<-	matrix(hpairs,1,2)
			} else {
				hpairs	<-	hpairs[order(hpairs[,1]),]
				}

				
		inc.slopes	<-	(hpairs[,2]-hpairs[,1])/(inc.points[hpairs[,2],2]-inc.points[hpairs[,1],2])
		
		inc.mle	<-	numeric()
		if(length(inc.slopes)>0){
			for (i in 1:length(inc.slopes)){
				inc.mle	<-	c(inc.mle, rep(inc.slopes[i], inc.points[hpairs[i,2],1]-inc.points[hpairs[i,1],1]))
				}
			}
		} else { inc.mle<-numeric()}


	mle		<- 	c(dec.mle, 0, inc.mle)

	H.vals 	<- 	cumsum(c(x[1],diff(x))*mle)
	S.vals 	<- 	exp(-H.vals)

 	mod.lik	<-    prod(dec.mle)*prod(inc.mle)
	lik		<-	prod(c(dec.mle, inc.mle)*S.vals[1:(n-1)])*S.vals[n]

	
	return(list(mle=mle, mod.lik=mod.lik, lik=lik, llh=log(lik)))
	}

## 	THIS FUNCTION FINDS THE MLE

find.mle	<- function(x){
	x		<-	sort(unique(x))
	n		<- 	length(x) #NEW
	llh	<-	rep(0,n)

	for (i in 1:n){
		llh[i]	<-	find.c.mle(x, i)$llh
		}
	
	loc		<-	which(llh==max(llh))

	if(length(loc)==1){
		mle	<- 	find.c.mle(x, loc)$mle
		}	else	{
			mle	<- numeric()
			}	

	unique	<-	FALSE
	if(length(loc)==1){
		unique	<- 	TRUE
		}		

	return(list(loc=loc, mle=mle, unique=unique))
	}