swh:1:snp:7d062e61dfe39b25a891852629fe3416f36554de
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
bjtest1d.R
bjtest1d <- function(y, d, x, beta) { ## depends on WKM( ), redistF( )
n <- length(y) ## dimension of x must be n x 1.
if ( length(x) != n ) stop("check dim of x")
if ( length(beta) != 1 ) stop("check dim of beta")
e <- y - beta * x
ordere <- order(e, -d)
esort <- e[ordere]
dsort <- d[ordere]
xsort <- x[ordere]
dsort[length(dsort)] <- 1 #last one as uncensored always
## use KM as F (need to be an n vector prob)
temp0 <- WKM(esort,dsort)
pold <- temp0$jump
temp <- redistF( y=esort, d=dsort, Fdist=pold )
weight <- temp$weight/n #the prob weight matrix
####pold <- colSums(weight) ##not needed, just let pold= WKM()$jump
A <- rep(0, n)
for (i in 1:n) if (dsort[i] == 1) {
A[i] <- sum( weight[1:i,i] * xsort[1:i] )/pold[i]
}
AA <- A[ dsort == 1 ]
myfun <- function(t, q) { t*q }
temp2 <- el.cen.EM(x=esort, d=dsort, fun=myfun, mu=0, q=AA)
pnew <- temp2$prob
logel1 <- temp0$logel
logel2 <- temp2$loglik
list(prob=pnew, logel=logel1, logel2=logel2, "-2LLR"=2*(logel1-logel2))
}