MeasureB <- function(B, lambda, ai, zi=numeric(0) ) {
#
# This function compute the measure beta, as the parameter of Diricrlet
# process prior. 
# ai is the partition, any given numbers, that need a measure computed.
# zi is the observed, uncensored data vector
#
     if(B<0) stop("weight B must >= 0")
     if(lambda<=0) stop("lambda must > 0")
     if(any(ai<=0)) stop("partition ai must > 0")

     times <- sort(ai)
     Malpha <- -B * diff(c(1, exp(-lambda * times), 0))

     if(length(zi)>0) {
     if(any(is.na(zi))) stop("NA's in the zi?")
     if(any(zi<0)) stop("negative zi?")
     Malpha <- Malpha-diff(c(length(zi),psum(t(outer(zi,times,FUN=">="))),0))
     }
     return(Malpha)
}

psum <- function(mat) {
         mat2 <- cbind(0, mat)
         as.vector(rowsum( t(mat2), c(0, rep(1, ncol(mat))))[2,])
}

Lemma1 <- function(B, lambda, ai, zi=numeric(0) ) {

BR <- B + length(zi)
N <- length(ai)-1 
fmu <- BR + (0:N)

Bi <- MeasureB(B, lambda, ai, zi)

fz <- (0:N) + cumsum(rev(Bi))[1:(N+1)]    

prod(fz/fmu) 
}