Some ideas of improvements:

To force the mean of the residual Kaplan-Meier = 0 
THis will also help define the last jump when the 
last obs is censored. (or may be this is auto matically 0 already?)




iter <- function(x, y, delta, beta)
{
# This function computes one iteration of the EM for 
# censored regression est. (Buckley-James estimator) It first order the
# data according to the residuals, then 
# compute conditional expectations, finally it compute a new beta by lm(). 
#
# Input: 
# x is a matrix of N rows (the covariates), 
# y is the censored response, a vector of length N.
# delta is a vector of length N.
# (delta =1 means (y) is uncensored. delta = 0 means
# the  (y) is right censored.)
# beta is the initial est. and is a vector of length = no. of column(x)
#
# Output: the new value of beta, 
#
        N <- length(delta)
        u <- x %*% beta
        res <- y - u               #### form the residuals
        niceorder <- order(res, - delta) # order the obs according to
        resorder <- res[niceorder]     # res, if tie then according to
        dorder <- delta[niceorder]  # delta value i.e. d=1 comes first
        dorder[N] <- 1              # added 2005 3/20
        uorder <- u[niceorder] 
        ystar <- y[niceorder]       # should I just let ystar <- delta ?
        xorder <- as.matrix(x[niceorder,])

temp <- WKM(x=resorder, d=dorder, zc=1:N)    # add  ( , zc=1:N )  2005/3

jifen <- rev( cumsum( rev(resorder * temp$jump)) )
Sresorder <- temp$surv

for (i in 1:N) if (dorder[i] == 0) {
           ystar[i] <- uorder[i] + jifen[i]/Sresorder[i]
}

return( lm( ystar ~ 0 + xorder )$coef )      ### the M-step
}





WKM <- function(x,d,zc=rep(1,length(d)),w=rep(1,length(d))) {
if (any((d != 0) & (d != 1)))
        stop("d must be 0(right-censored) or 1(uncensored)")
temp <- Wdataclean3(x,d,zc,w)
dd <- temp$dd
ww <- temp$weight
dd[length(dd)] <- 1

allrisk <- rev(cumsum(rev(ww)))
survP <- cumprod( 1 -  (dd*ww)/allrisk )
jumps <- -diff( c(1, survP) )

logel <- sum(ww[dd==1]*log(jumps[dd==1])) + sum(ww[dd==0]*log(survP[dd==0]))

list(times=temp$value, jump=jumps, surv=survP, logel=logel)

}
################################################################
##This function compute a weighted Kaplan-Meier estimator
## x = times, d = censoring status, w = weights 
################################################################