# R code to find the MLE of a decreasing mass function on {0,1,2, ...}.  The input of the function  is the vector 
# $\widehat p_n$ and the output is $\gren(\widehat p_n)$.  If \texttt{plot=TRUE}, then a visual representation of 
# the LCM is also given.

find.gren  <-  function(p, plot=FALSE){
    n      <-  length(p)
    pts    <-  cbind(c(0:n,n),c(0,cumsum(p),0))

    hpts   <-  chull(pts)
    hpts   <-  c(hpts, hpts[1])
	
    if(plot==TRUE){
    plot(pts, cex = 1, pch=19, ylab="", xlab="", xlim=c(0,n), ylim=c(0, 1))
    lines(pts[hpts, ])
    points(pts[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(pts[,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]),]
                }

    s.num   <-  pts[hpairs[,1],2]-pts[hpairs[,2],2]
    s.denom <-  pts[hpairs[,1],1]-pts[hpairs[,2],1]				
    slopes  <-  s.num/s.denom
    h.hat   <-  numeric()
    for (i in 1:length(slopes)){
        h.hat	<-	c(h.hat, rep(slopes[i], hpairs[i,1]-hpairs[i,2]))
        }

    return(h.hat)
    }