https://github.com/cran/emplik
Tip revision: 6d2aaee84bc7d13a97cf480fbaf9e89663bcb3cd authored by Mai Zhou on 09 January 2009, 00:00:00 UTC
version 0.9-5
version 0.9-5
Tip revision: 6d2aaee
el.cen.EM.R
el.cen.EM <- function(x,d,fun=function(t){t},mu,maxit=25,error=1e-9, ...) {
xvec <- as.vector(x)
nn <- length(xvec)
if(nn <= 1) stop ("Need more observations")
if(length(d)!=nn) stop("length of x and d must agree")
if(any((d!=0)&(d!=1)&(d!=2)))
stop("d must be 0(right-censored) or 1(uncensored) or 2(left-censored)")
if(!is.numeric(xvec)) stop("x must be numeric")
if(length(mu)!=1) stop("check the dim of constraint mu")
temp <- Wdataclean2(xvec,d)
x <- temp$value
d <- temp$dd
w <- temp$weight
###### change the last obs. among d=1,0, so that d=1
###### change the first obs. among d=1,2 so that d=1
###### this ensures we got a proper dist. for NPMLE. (no mass escape)
INDEX10 <- which(d != 2)
d[ INDEX10[length(INDEX10)] ] <- 1
INDEX12 <- which(d != 0)
d[ INDEX12[1] ] <- 1
xd1 <- x[d==1]
if(length(xd1) <= 1) stop("need more distinct uncensored obs.")
funxd1 <- fun(xd1, ...)
xd0 <- x[d==0]
xd2 <- x[d==2]
wd1 <- w[d==1]
wd0 <- w[d==0]
wd2 <- w[d==2]
m <- length(xd0)
mleft <- length(xd2)
##############################################
#### do the computation in 4 different cases.#
##############################################
if( (m>0) && (mleft>0) ) {
pnew <- el.test.wt(funxd1, wt=wd1, mu)$prob
n <- length(pnew)
k <- rep(NA, m)
for(i in 1:m) { k[i] <- 1+n - sum( xd1 > xd0[i] ) }
kk <- rep(NA, mleft)
for(j in 1:mleft) { kk[j] <- sum( xd1 < xd2[j] ) }
num <- 1
while(num < maxit) {
wd1new <- wd1
sur <- rev(cumsum(rev(pnew)))
cdf <- 1 - c(sur[-1],0)
for(i in 1:m)
{wd1new[k[i]:n] <- wd1new[k[i]:n] + wd0[i]*pnew[k[i]:n]/sur[k[i]]}
for(j in 1:mleft)
{wd1new[1:kk[j]] <- wd1new[1:kk[j]] + wd2[j]*pnew[1:kk[j]]/cdf[kk[j]]}
temp8 <- el.test.wt(funxd1, wt=wd1new, mu)
pnew <- temp8$prob
lam <- temp8$lam
num <- num +1
}
logel <- sum(wd1*log(pnew)) + sum(wd0*log(sur[k])) + sum(wd2*log(cdf[kk]))
tempDB <- WCY(x=x,d=d,wt=w)
logel00 <- tempDB$logEL
funNPMLE <- sum( fun(tempDB$time)* tempDB$jump )
}
if( (m>0) && (mleft==0) ) {
temp3 <- WKM(x=x,d=d,w=w)
logel00 <- temp3$logel
funNPMLE <- sum(funxd1 * temp3$jump[temp3$jump > 0])
# now the iteration
pnew <- el.test.wt(x=funxd1, wt=wd1, mu=mu)$prob
n <- length(pnew)
k <- rep(NA,m)
for(i in 1:m) { k[i] <- 1 + n - sum( xd1 > xd0[i] ) }
num <- 1
while(num < maxit) {
wd1new <- wd1
sur <- rev(cumsum(rev(pnew)))
for(i in 1:m)
{wd1new[k[i]:n] <- wd1new[k[i]:n] + wd0[i]*pnew[k[i]:n]/sur[k[i]]}
temp9 <- el.test.wt(funxd1, wt=wd1new, mu)
pnew <- temp9$prob
lam <- temp9$lam
num <- num +1
}
sur <- rev(cumsum(rev(pnew)))
logel <- sum( wd1*log(pnew)) + sum( wd0*log(sur[k]) )
}
if( (m==0) && (mleft>0) ) {
kk <- rep(NA, mleft)
for(j in 1:mleft) { kk[j] <- sum( xd1 < xd2[j] ) }
pnew <- el.test.wt(funxd1, wt=wd1, mu)$prob
n <- length(pnew)
num <- 1
while(num < maxit) {
wd1new <- wd1
cdf <- cumsum(pnew)
for(j in 1:mleft)
{wd1new[1:kk[j]] <- wd1new[1:kk[j]] + wd2[j]*pnew[1:kk[j]]/cdf[kk[j]]}
temp7 <- el.test.wt(funxd1, wt=wd1new, mu)
pnew <- temp7$prob
lam <- temp7$lam
num <- num +1
}
logel <- sum( wd1*log(pnew)) + sum( wd2*log( cdf[kk] ) )
dleft <- d
dleft[dleft==2] <- 0
templeft <- WKM(x= - x, d=dleft, w=w) ### bug fix 3/2008
logel00 <- templeft$logel ### substitute a left WKM() ???
funNPMLE <- NA
}
if( (m==0) && (mleft==0) ) {
funNPMLE <- sum( funxd1 * wd1/sum(wd1) )
logel00 <- sum( wd1*log( wd1/sum(wd1) ) )
temp6 <- el.test.wt(funxd1, wt=wd1, mu)
pnew <- temp6$prob
lam <- temp6$lam
logel <- sum( wd1*log(pnew) )
}
# get ready for exit
tval <- 2*(logel00 - logel)
list(loglik=logel, times=xd1, prob=pnew, funMLE=funNPMLE, lam=lam,
"-2LLR"=tval, Pval= 1-pchisq(tval, df=1) )
}