.packageName <- "NPBayes"
Lemma1 <- function(B, theta, ai, zi=numeric(0)) {
          if(B<=0) stop("weight B must be > 0")
          if(theta<=0) stop("theta must be > 0")
          if(any(ai<0)) stop("partition ai must be >= 0")
          if(all(ai==0)) return(1)
             else { ai <- ai[ai>0]
                    BR <- B + length(zi)
                    N <- length(ai)-1
                    fmu <- BR + (0:N)
                    Bi <- MeasureB(B, theta, ai, zi)
                    fz <- (0:N) + cumsum(rev(Bi))[1:(N+1)]
                    return(prod(fz/fmu))
                  }
}

### file NPBayes/R/npbayes.R
### copyright (C) 2004-6 Mai Zhou

Allends <- function(lefts, rights) {
#  lefts   are the left end points of intervals,
#  rights  are the right end points of intervals.
if(length(lefts)!=length(rights)) 
                   stop("Lefts and Rights must have same length")
if(any(lefts>=rights)) stop("Rights must > Lefts")
if(length(lefts)>0)
 return(t(1-binary4(length(lefts)))*lefts+t(binary4(length(rights)))*rights)
}

Allsigns <- function(k) {
# generate the signs for the column of Allends().
x1<-c(1,-1)
x2<-c(x1, -x1)
x3<-c(x2, -x2)
x4<-c(x3, -x3)
x5<-c(x4, -x4)
if(k==0) return(1)
if(k==1) return(x1)
if(k==2) return(x2)
if(k==3) return(x3)
if(k==4) return(x4)
if(k==5) return(x5)
if(k>5) {xold<-x5
         num<-5
         while(k>num) {num<-num+1
                       xnew<-c(xold,-xold)
                       xold<-xnew
                      }
         return(xnew)
        }
}

MeasureB <- function(B, theta, ai, zi=numeric(0) ) {
# This function computes the measure beta, as the parameter for Dirichlet
# process prior.
# ai is the partition, any given >0 numbers, that need a measure computed.
# zi is the observed, uncensored data vector.
#
     times <- sort(ai)
     Malpha <- -B * diff(c(1, exp(- theta * 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),rowSums(t(outer(zi,times,FUN=">="))),0))
     }
     return(Malpha)
}

binary4 <- function(k){ as.matrix( rev(expand.grid(rep(list(0:1),k)))) }

### file NPBayes/R/npbayes.R
### copyright (C) 2004-6 Mai Zhou

NPBayes <- function(B, theta, u, uncen=numeric(0), rightcen=numeric(0),
                                   lefts=numeric(0), rights=numeric(0)) {
# B and theta are required inputs, they are the parameters for the
# Dirichlet process prior. i.e. the measure \alpha (t, \infty) =
# B \exp( -\theta t). 
# u (required), a scaler, is the time we want to have 1-F(u) computed.
# uncen, rightcen are the uncensored and right censored observations.
# lefts, rights are the end points of interval censored observations.
# If there are no uncensored observations we should say uncen=numeric(0) etc.
# The left censored observations can be handled by treating them as
# interval censored: [0, a) == left censored at a.
#
# Due to rounding errors and loss of significant digits, if there are a lot
# of interval censored data, the result can be untrustworthy.

  if(length(u)>1) stop("currently u can only be scalar")
  matends <- Allends(lefts=lefts, rights=rights)
  dimnames(matends) <- NULL
  vecsigns <- Allsigns(length(lefts))
  M <- length(vecsigns)
  numer <- denom <- rep(NA, M)

 for(i in 1:M) {
    partit <- c(rightcen, matends[,i])
    numer[i] <- Lemma1(B=B, theta=theta, ai=c(u,partit), zi=uncen)
    denom[i] <- Lemma1(B=B, theta=theta, ai=partit, zi=uncen)
    }
  return(sum(numer*vecsigns)/sum(denom*vecsigns))
}

# To do list: (1) To make the choice of alpha measure more flexible.
# Instead of the type B exp(-theta t), let the user supply their
# own measure. (2) To remove the restriction of non-negative data.
# so the alpha measure can be on the whole line instead of positive line.
### NPBayes for truncated and censored data
### GPL copyright (C) 2004-6 Mai Zhou

NPBayesT <- function(B, theta, u, y=numeric(0), x=numeric(0), 
                                           status=numeric(0) ) {
# B and theta are required inputs, they are the parameters for the
# Dirichlet process prior. i.e. the measure \alpha (t, \infty) =
# B \exp( -\theta t). 
# u (required), a scaler, is the time we want to have the Bayes est.
#                    1- \hat F(u) computed.
# y, is the vector of truncation times. 
# x, is the observed (may be right censored) lifetimes .
# status, is the indicator of right censoring.
#
  if(length(u)>1) stop("currently u can only be a scalar")
  if(B <= 0) stop("B must be > 0")
  if(theta <= 0) stop("theta must be > 0")

  fenzi <- NewLemma1(B=B, theta=theta, u=u, y=y, x=x, d=status)
  fenmu <- NewLemma1(B=B, theta=theta, y=y, x=x, d=status)
  return(fenzi/fenmu) 
}

# To do list: (1) same as in NPBayes (2) to save calculation by
# omitting those term that will be canceled in the ratio of fenzi/fenmu.
### copyright (C) 2004 Mai Zhou 
# Compute an expectation of a likelihood. 
# B and theta are required inputs, they are the parameters for the
# Dirichlet process prior. i.e. the measure \alpha (t, \infty) =
# B \exp( -\theta t). 
# y is the (left) truncation times.
# x is the (possibly) right censored observations.
# d is the censoring status.
#   if there is a u, then the likelihood include a term P([u,\infty)).

NewLemma1 <- function(B, theta, u = numeric(0), y, x, d) {
if(B<=0) stop("weight B must be > 0")
if(theta<=0) stop("theta must be > 0")
n <- length(y)
if(length(x) != n) stop("length of x and y must agree")
if(length(d) != n) stop("length of d and y must agree")
if(any((d!=0)&(d!=1))) stop("d must be 0/1's for censor/not-censor")
if( any(x <= y) ) stop("x must > y, as y is the entry time")

if(length(u) == 0) 
mydata <- matrix(c(y,x,rep(-1,n),rep(1,n),rep(9,n),d), 
                            nrow=3, ncol=2*n, byrow=TRUE)

if(length(u) == 1) 
mydata <- matrix(c(u,y,x, 1,rep(-1,n),rep(1,n), 0,rep(9,n),d), 
                            nrow=3, ncol=(2*n +1), byrow=TRUE)

if(length(u) > 1) stop("u can only be a scalar")
zorder <- order(mydata[1, ])
sortdata <- mydata[ ,zorder]
## first row is the times y and x. Second row is the power in Likelihood.
## Third row is the censoring indicator, 0=right censor, 1=not censor.
if(all(sortdata[1,]==0)) return(1) 
if(any(sortdata[1,] < 0)) stop("y,x must > 0")
else { ns <- rev(cumsum(rev(sortdata[2,])))
       ns <- c(ns[-1], 0) 
       ms <- as.numeric(sortdata[2,]==-1)
       Bi <- B * exp(- theta * sortdata[1,]) 
       fz <- ( Bi + ns - ms )^sortdata[2,]
       fz <- fz[sortdata[3,] != 1] 
       fz <- prod(fz)
       if(length(u) == 1) fz <- fz/B
       return(fz)
      }
}

