https://github.com/cran/emplik
Tip revision: 1a00a3dce7a5ac3854ad4978598dbc8ab597c49f authored by Mai Zhou on 23 October 2010, 00:00:00 UTC
version 0.9-6
version 0.9-6
Tip revision: 1a00a3d
el.cen.EM2.R
el.cen.EM2 <- function(x,d,xc=1:length(x),fun,mu,maxit=25,error=1e-9,...){
####
#### xc is collaps control: if for index i and j, (x[] d[]) are identical
#### should they be merged into one obs. with weight 2?
#### if xc[i] != xc[j] then they do not merge.
#### For regression applications, you should not merge -- use
#### xc = 1:length(x) .
#### The order of censoring (re-distribute weights)
#### is based on x, not fun(x).
####
xvec <- as.vector(x)
d <- as.vector(d)
mu <- as.vector(mu)
xc <- as.vector(xc)
n <- length(d)
if (length(xvec)!=n) stop ("length of d and x must agree")
if (length(xc)!=n) stop ("length of xc and d must agree")
if(n <= 2*length(mu)+1) stop ("Need more observations")
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(!is.numeric(mu)) stop("mu must be numeric")
funx <- as.matrix(fun(xvec, ...)) ##get the matrix before sorting xvec
pp <- ncol(funx)
if(length(mu)!=pp) stop("length of mu and ncol of fun(x) must agree")
temp <- Wdataclean5(z=xvec,d,zc=xc, xmat=funx) ## collaps control zc
x <- temp$value
d <- temp$dd
w <- temp$weight
funx <- temp$xxmat ## make sure funx and x are in same order
###### 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 CDF for NPMLE. (no mass escapes)
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 <- funx[d==1,]
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.wt2(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] ) }
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.wt2(x=funxd1, wt=wd1new, mu=mu)
pnew <- temp8$prob
lam <- temp8$lambda
num <- num +1
}
logel <- sum(wd1*log(pnew)) + sum(wd0*log(sur[k])) + sum(wd2*log(cdf[kk]))
logel00 <- WCY(x=x,d=d,wt=w)$logEL #### should I add zc?
}
if( (m>0) && (mleft==0) ) {
pnew <- el.test.wt2(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.wt2(x=funxd1, wt=wd1new, mu=mu)
pnew <- temp9$prob
lam <- temp9$lambda
num <- num+1
}
sur <- rev(cumsum(rev(pnew)))
logel <- sum( wd1*log(pnew)) + sum( wd0*log(sur[k]) )
logel00 <- WKM(x=x,d=d,zc=xc,w=w)$logel
}
if( (m==0) && (mleft>0) ) {
kk <- rep(NA, mleft)
for(j in 1:mleft) { kk[j] <- sum( xd1 < xd2[j] ) }
pnew <- el.test.wt2(x=funxd1, wt=wd1, mu=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.wt2(x=funxd1, wt=wd1new, mu=mu)
pnew <- temp7$prob
lam <- temp7$lambda
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, zc=xc, w=w) ### bug fix 3/2008
logel00 <- templeft$logel ### ?? do I need a left WKM ??
}
if( (m==0) && (mleft==0) ) {
num <- 0
temp6 <- el.test.wt2(x=funxd1, wt=wd1, mu)
pnew <- temp6$prob
lam <- temp6$lambda
logel <- sum( wd1*log(pnew) )
logel00 <- sum( wd1*log( wd1/sum(wd1) ) )
}
tval <- 2*(logel00 - logel)
list(loglik=logel, times=xd1, prob=pnew, lam = lam, iters=num,
"-2LLR"=tval, Pval= 1-pchisq(tval, df=length(mu)) )
}