%From chen@ms.uky.edu  Mon Aug 28 15:53:13 2000
%Dr. Zhou, %Here is my program for right censored data
%using your new idea.  Kun

#qprcinv <- function(x,d,w0=1/length(d[d==1]),mu,error=1e-8) {
#    
#   #**********************************************************   
#   # This program computes the empirical loglikehood function
#   # for right censored data (x,d) with constraint: 
#   # sum(w*x*d) = mu ,  sum(w*d)=1  and  I(d=0)*wj = sum (wi*d)
#   #
#   # Required Inputs:
#   #     x   right-censored survival times. a vector
#   #     d   censoring indicator vector 1-uncensor, 0-censor
#   #     mu  constrained mean value of the survival distribution, scalar 
#   # Optional inputs:
#   #     w0     initial weight value, a vector. 
#   #     error  iteration stopping control standard, positive scalar.
#
#   # Output
#   # loglikfn: loglikelihood function
#   # weigth  : weight under constraint
#   # m       : the number of iteration
#   #**********************************************************
# 
#   # Before computation, we should clean data first. 
#   # If the last obs. is censored, we should change it to uncensored
#   # for easy calculation.
#
#To: mai Zhou <mai@ms.uky.edu>
#From: Kun Chen <chen@ms.uky.edu>
#Subject: 
#Mime-Version: 1.0
#Content-Type: multipart/mixed;
#Status: OR
#
#Hi, Dr. Zhou,
#
#I think I finally finished the right censored programs.
#I tested on the data from Thomas & Grunkemier.
#It works.
#
#Here are these program package.
#
#Kun 
#
#fun <- function(x,t){as.numeric(x<=t)}

qprcinvfun <- function(x,d,w0=NA,mu,fun,t,
              identical=rep(0,length(x)), error=1e-2)
{
   
 #**********************************************************   
   # This program compute the loglikehood function
   # for right censored data with constraint: 
   # sum(w*fun(x)) = mu

   # Input
   # x  <- right-censored life time data vector
   # d  <- censoring indicator vector
   # mu <- constrained value
   # w0 <- initial weight value
   # error <- iteration stopping control standard

   # Output
   # loglikfn: loglikelihood function
   # weigth  : weight under constraint
   # m       : the number of iteration

 #**********************************************************
   n <- length(x)
   if(w0=="NA")
   {
      niceorder <- order(x,-d)
      sortx <- x[niceorder]
      sortd <- d[niceorder]
      
   #Change the last right censored obs as uncensored.
   
      last <- length(sortd)
      if(sortd[last] == 0) 
      {
         i <- last
         while(sortd[i] != 1 && i > 0) 
         {
            if(sortd[i] == 0 && sortx[i] == sortx[last])
            {               
              sortd[i] <- 1
            }
            i <- i - 1
         }
           
      }
   
      source("u:/research/QP/testmean/KM.r")
      KM <- KapMei(sortx,sortd)
      sortw0 <- KM$jump
  #    print("KM")
  #    print(KM)
     
  
   }
   else
   {
   # Before computation, we should clean data first.
   # If the last obs. is censored, we should change to uncensored
   # for easy calculation.

   
      source("u:/research/QP/testmean/cleanrc.txt")
      data <- cleanrc(x,d,n,w0,identical)
      sortx <- data$sortx
      sortd <- data$sortd
      sortw0 <- data$sortw0
      weight <- data$weight
    
      brtie <- c(-9999,diff(sortx))
      sortx[brtie==0] <- sortx[brtie==0]+0.00001*c(1:length(brtie[brtie==0]))
   }

   # If all obs. are uncensored, we use another program to
   # calculate the loglikelihood.
   if(all(sortd==1))
   {
    
     source("u:/research/QP/testmean/qpuncen.txt")
     qpuncen(sortx,sortw0,mu,n,fun,error)
     
   }
   else
   {
  
   uncensortx <- sortx[sortd==1]
 
   uncensortw0 <- sortw0[sortd==1]
 #  print("uncensortw0")
 #  print(uncensortw0)
 #  print("uncensortx")
 #  print(uncensortx)
 
   k <- length(uncensortx)
   xbar <- sum(fun(uncensortx,t)*uncensortw0)  
  
  
 #************* begin initial calculation***********************

   # get vector dvec which is the first derivative vector
  
   dvec01 <- weight[sortd==1]/uncensortw0  
   rk <- rank(sortx)
   uncenrk <- rk[sortd==1]
   cenrk <- rk[sortd==0]
   mm <- length(cenrk)
   dvec02 <- rep(0,mm)
#   print("uncenrk")
#   print(uncenrk)
#   print("cenrk")
#   print(cenrk)
  
  
   len <- length(sortw0)
   for(j in 1:mm)
   {   
      vv <- sum(sortw0[cenrk[j]:len]) 
      sortw0[cenrk[j]] <- vv
      dvec02[j] <- weight[cenrk[j]]/vv   
   }   
 
   dvec0 <- rep(0,len)
   dvec0[sortd==1] <- dvec01
   dvec0[sortd==0] <- dvec02
#   print("dvec0")
#   print(dvec0)

  # get matix Dmat which is 2nd derivative matrix.
  
  Dmat0 <- diag(sqrt(weight)/dvec0)
#    print("Dmat0")
#  print(Dmat0)
 
   mat <- matrix(rep(sortd,mm),ncol=mm)
   for(i in 1:mm)
   {
      mat[,i][1:cenrk[i]] <- 0
      mat[,i][cenrk[i]] <- -1
   }
 
  # get constraint matrix Amat
   Amat <- as.matrix(cbind(sortd,fun(sortx,t)*sortd, mat))
#   print("Amat")
#   print(Amat)
  
  # get constraint vector bvec 
   bvec0 <- c(0,mu-xbar,rep(0,mm))
#   print("bvec0")
#   print(bvec0)
  
  # Use solve.QP to maximize the loglikelihood function
   value0 <- solve.QP(Dmat0,dvec0,Amat,bvec0,meq=mm+2,factorized=T)
   
  w <- sortw0+value0$solution
 
   if(any(w<=0)) 
       stop("Weight should be positive value. There is no probability satisfying this constraint")

 #**************end initial calculation *******************************************                       

 #**************begin iteration ***************************************************
   
   # update vector Dmat, dvec and bvec after initial calculation
   bvec <- rep(0,mm+2)
   diff <- 10

   m <- 0
   while(diff > error)
   {
     dvec <- weight/w
  
   
     # get matix Dmat
     Dmat <- diag(w/sqrt(weight)) 
    
     value0 <- solve.QP(Dmat,dvec,Amat,bvec,meq=mm+2, factorized=T)
   
     w <- w + value0$solution

     diff <- sum(value0$solution^2)
     m <- m+1
   }
   loglikfn <- sum(log(w))
  
   list(loglikfn=loglikfn, weight=w, iteration=m, error=error)

 }
   
}

--=====================_11560913==_
Content-Type: text/plain; charset="us-ascii"
Content-Disposition: attachment; filename="cleanrc.txt"

cleanrc <- function(x,d,n,w0,identical=rep(0,length(x)))
{
   if(n <= 2) stop ("Need more observation")
   if(length(d) != n) stop ("length of x and d must agree")
   if(length(w0)!=n)  stop("length of x and w0 must agree")
   if(any((d!=0)&(d!=1)))
      stop ("d must be 0(right-censored), 1(uncensored)")
   if(all(d==0))
      stop ("All Obs. are censored. Invalid Data")
   if(!is.numeric(x))
     stop("x must be numeric values --- observed times")   
   if(sum(w0)<0.999)
     stop("Invalid Initial Weight")
   if(any(w0<0))
     stop("Weight must be nonegative")

  
   niceorder <- order(x,  - d)   # order them as z increase. When z 
   sortx <- x[niceorder]         # has ties, use d=2, 1, 0 to order 
   sortd <- d[niceorder]         # the tied z values. 
   sortw0 <- w0[niceorder]
   dupsortx <- duplicated(sortx)  #see if there is dup in z. But even
   argdiff <- c(1, diff(sortd))  #z's tie, if d is diff, it's not a tie
   dupsortx[argdiff != 0] <- F    #seems I need not dupsortz ==T &
   dupsortx[identical != 0] <- F  #also, do not collaps if identical !=0

   sortx <- sortx[dupsortx != T]   # get the unique values of sortz and 
   sortd <- sortd[dupsortx != T]   # sortd. the weight or duplicatility
   sortw0 <- sortw0[dupsortx != T] # w[i] will need to be computed later in C function

   count <- (1:length(dupsortx))[dupsortx != T]
   weight <- diff( c(count, (1+length(dupsortx)) ) )

   #Change the last right censored obs as uncensored.
   
   last <- length(sortd)
   if(sortd[last] == 0) 
   {
       i <- last
       while(sortd[i] != 1 && i > 0) 
       {
          if(sortd[i] == 0 && sortx[i] == sortx[last])
          {               
              sortd[i] <- 1
             
          }
          i <- i - 1
       }
           
   }
   
   sortw0[sortd==1] <- weight[sortd==1]/sum(weight[sortd==1])
   list(sortx=sortx, sortd=sortd, sortw0=sortw0,weight=weight)

}
--=====================_11560913==_
Content-Type: text/plain; charset="us-ascii"
Content-Disposition: attachment; filename="QPuncen.txt"

fun <- function(x,t){as.numeric(x<=t)}

qpuncen <- function(x, w0, mu,fun,t,error=1e-8)
{   
   n <- length(x)
   niceorder <- order(x)
   sortx <- x[niceorder]
   sortw0 <- w0[niceorder]   
   xbar <- sum(fun(sortx,t)*sortw0)    
   print(xbar) 
   mat0a <- diag(sortw0)   
   mat0b <- 1/sortw0   
   mat0c <- cbind(rep(1,n),fun(sortx,t))   
   print("mat0c")
   print(mat0c)
   mat0d <- c(0, mu-xbar)   
   print("mat0d")
   print(mat0d)
   value0 <- solve.QP(mat0a,mat0b, mat0c, mat0d, meq=2,factorized=T)  
   w <- sortw0 + value0$solution 

   #iteration  
   matd <- c(0,0)  
   diff <- 10  
   m <- 0  
   while(diff>error)  
   {    
      mata <- diag(w)    
      matb <- 1/w    
      value <- solve.QP(mata,matb,mat0c,matd,meq=2,factorized=T)    
      w <- w+value$solution    
      diff <- sum(value$solution^2)    
      m <- m+1 
   }    
   print(sum(w))
   print(sum(w[14:21]))   
   likefn <- sum(log(n*w))     
   list(likefn=likefn, Weight=w, iteration=m, error=error)   
 
}