samplemaker_gamma <- function(n=100,shape=2,scale=1,end=4 ){
x <- rgamma(n,shape = 2,scale = 1)
t <- runif(n,0,4)
delta <- 0.5
M <- 1
observedx <- x[which(x>t)]
observedt <- t[which(x>t)]
data <- data.frame(observedx,observedt)
status <- ifelse( observedx<=observedt+M*delta, 3  , 0) #maybe change
qi <- ifelse(status==3,  observedt  ,observedt+M*delta   )
pi <- ifelse( status==3,  observedt+M*delta ,  Inf                 )
datause <- data.frame(observedx,observedt ,Li=qi,Ui=pi,status,entry=observedt)
return(datause)
}
 
cria.tau <- function(data){
l <- data$Li
r <- data$Ui
tau <- sort(unique(c(l,r[is.finite(r)])))
#tau <- c( 0,tau,100000)
return(tau)
}
makeinterval <- function(data){
x <- matrix( c(0,rep(cria.tau( data ),rep(2,length( cria.tau( data )  ))) ,100000 )
,ncol = 2 ,byrow = T)
q <- x[,1]
p <- x[,2]
data.frame(q,p)
}
 
make_A_turnbull <- function(datacheck){
  qj<- makeinterval(datacheck)$q
  pj<- makeinterval(datacheck)$p
  li <- datacheck$Li
  ri <- datacheck$Ui
  n <- length(li)
  m <- length(pj)
  A <- matrix(0,nrow = n,ncol = m)
  intmapL <- signif(qj,4)
  intmapR <- signif(pj,4)
  k <- dim(A)[[2]]
  Lbracket <- rep("(", k)
  Rbracket <- rep(")", k)
  intname <- paste(Lbracket, intmapL, ",", intmapR, Rbracket,
  sep = "")
 
  intmapL2 <- signif(datacheck$Li, 4)
  intmapR2 <- signif(datacheck$Ui, 4)
  k2 <- dim(A)[[1]]
  Lbracket2 <- rep("(", k2)
  Rbracket2 <- rep(")", k2)
  intname2 <- paste(Lbracket2, intmapL2, ",", intmapR2, Rbracket2,
  sep = "")
  for(j in 1:m){
    for( i in 1:n){
      if( li[i]<=qj[j]&&ri[i]>=pj[j]      ){
        A[i,j] <- 1
      }
      else{A[i,j] <- 0}
    }
  }
 
colnames(A) <- intname
rownames(A) <- intname2
return(A)
}
 
前面的東西不用管 只是製作A矩陣必要的CODE

set.seed(5)
datatest <- samplemaker_gamma(n=7)
print("datatest")
datatest
 
A <- make_A_turnbull(datatest)
### in this case you can see their are some zero in the matrix
A
 
n<-nrow(A)     #number of observation
m<-ncol(A)     #number of interval
s <- rep(   1/dim(A)[2],  dim(A)[2] )
###  my code to compute expected death
 
den <- rep(0,n)
densum <- rep(0,n)
d<- rep(0,m)
 
for(j in 1:m){  
 
for(i in 1:n){
 
  if(A[i,j]>0){ ##因為有0在A矩陣的元素中 所以我想到這樣做 但有更好的做法嗎
 
  den[i]<- sum(A[i,]*s)
  densum[i]<- A[i,j]*s[j]/den[i]
 
  }
}
d[j]<- sum( densum )
}
d