https://github.com/cran/emplik
Tip revision: e2f16b0adc0f6df124c7ee126345f66d8d2c7961 authored by Mai Zhou on 09 October 2011, 00:00:00 UTC
version 0.9-7
version 0.9-7
Tip revision: e2f16b0
el.ltrc.EM.R
el.ltrc.EM <- function(y,x,d,fun=function(t){t},mu,maxit=30,error=1e-9) {
yvec <- as.vector(y)
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)))
stop("d must be 0(right-censored) or 1(uncensored)")
if(!is.numeric(xvec)) stop("x must be numeric")
if(length(mu)!=1) stop("check the dim of constraint mu")
yvec <- yvec[yvec > -Inf]
N <- length(yvec)
if ( N == 0 ) {
temp1 <- el.cen.EM(xvec,d,fun=fun,mu=mu)
WILKS <- temp1$"-2LLR"
pnew <- temp1$prob
}
temp <- Wdataclean2(xvec,d)
x <- temp$value
d <- temp$dd
w <- temp$weight
###### change the last obs.'s censoring indicator, so that d=1
###### this ensures we got a proper df for NPMLE. (no mass escape)
d[length(d)] <- 1
xd1 <- x[d==1]
funxd1 <- fun(xd1)
n <- length(xd1)
if(n <= 1) stop("need more distinct uncensored obs.")
xd0 <- x[d==0]
wd1 <- w[d==1]
wd0 <- w[d==0]
mright <- length(xd0)
if ( mright == 0 ) {
temp2 <- el.trun.test(yvec,xd1,fun=fun,mu=mu)
WILKS <- temp2$"-2LLR"
pnew <- temp2$NPMLEmu
}
if( mright>0 ) {
p0 <- LTRC(x,d,w,yvec)$survjump
pnew <- p0
k <- rep(NA,mright)
for(i in 1:mright) { k[i] <- 1 + n - sum( xd1 > xd0[i] ) }
indi <- function(u,v){ as.numeric(u > v) }
uij <- outer(xd1,yvec,FUN="indi")
num <- 1
while(num <= maxit) {
wd1new <- wd1
##########right censor
sur <- rev(cumsum(rev(pnew)))
for(i in 1:mright)
{wd1new[k[i]:n] <- wd1new[k[i]:n] + wd0[i]*pnew[k[i]:n]/sur[k[i]]}
#########left truncated
pvec <- as.vector( pnew %*% uij )
wd1new <- wd1new + as.vector(rowsum(t(pnew*(1-uij))/pvec, group=rep(1,N)))
######### weighted computation
pnew <- el.test.wt(funxd1, wt=wd1new, mu=mu)$prob
num <- num +1
}
sur <- rev(cumsum(rev(pnew)))
pvec <- as.vector( pnew %*% uij )
logel <- sum(wd1*log(pnew))+sum(wd0*log(sur[k]))-sum(log(pvec))
sur0 <- rev(cumsum(rev(p0)))
pvec0 <- as.vector( p0 %*% uij )
logel0 <- sum(wd1*log(p0))+sum(wd0*log(sur0[k]))-sum(log(pvec0))
WILKS <- 2*(logel0 - logel)
}
list(times=xd1, prob=pnew, "-2LLR"= WILKS, Pval=1-pchisq(WILKS,df=1))
}