New.el2.cen.EMs <- function (x, y, fun = function(x, y) {x >= y}, mean = 0.5, maxit = 25) 
{
    xvec <- as.vector(x)
    yvec <- as.vector(y)
    nx <- length(xvec)
    ny <- length(yvec)
    if (nx <= 1) 
        stop("need more observations in x")
    if (ny <= 1) 
        stop("need more observations in y") 
    if (!is.numeric(xvec)) 
        stop("x must be numeric")
    if (!is.numeric(yvec)) 
        stop("y must be numeric")
    if (length(mean) != 1) 
        stop("mean must have dimension 1")
    temp1x <- myWdataclean2(z=xvec, d=rep(1,nx))
    temp1y <- myWdataclean2(z=yvec, d=rep(1,ny))
    x <- temp1x$value
    y <- temp1y$value
    dx <- temp1x$dd
    dy <- temp1y$dd
    wx <- temp1x$weight
    wy <- temp1y$weight
    xindex10 <- which(dx != 2)
    yindex10 <- which(dy != 2)
    dx[xindex10[length(xindex10)]] <- 1
    dy[yindex10[length(yindex10)]] <- 1
    xindex12 <- which(dx != 0)
    yindex12 <- which(dy != 0)
    dx[xindex12[1]] <- 1
    dy[yindex12[1]] <- 1
    
    xd1 <- x[dx == 1]
    wxd1 <- wx[dx == 1]
    yd1 <- y[dy == 1]
    wyd1 <- wy[dy == 1]
    
    if (length(xd1) <= 1) stop("need more distinct uncensored x obs.")
    if (length(yd1) <= 1) stop("need more distinct uncensored y obs.")
    
    nx1 <- length(xd1)
    ny1 <- length(yd1)   
    nx <- length(x)
    ny <- length(y)   
   

    logelx00 <- sum(wxd1 * log(wxd1/sum(wxd1)))
    jumpxu <- wxd1/sum(wxd1)    
   
    logely00 <- sum(wyd1 * log(wyd1/sum(wyd1)))
    jumpyu <- wyd1/sum(wyd1)
   
    logel00 <- logelx00 + logely00
    indic <- matrix(NA, nrow = nx1, ncol = ny1)
    for (i in 1:nx1) {
        for (j in 1:ny1) {
            indic[i, j] <- fun(xd1[i], yd1[j])
        }
    }
    indicmat <- indic - mean
    funNPMLE = as.vector(jumpxu %*% indic %*% jumpyu)
    muvec <- jumpxu
    nuvec <- jumpyu
    

    num <- 1
    logvec <- rep(0, maxit)
    while (num <= maxit) {
        temp3 <- New.el2.test.wts(wxd1, wyd1, muvec, nuvec, indicmat)
        muvec <- temp3$jumpu
        nuvec <- temp3$jumpv
        logelx <- sum(wxd1 * log(muvec))
        logely <- sum(wyd1 * log(nuvec))
        logel <- logelx + logely
        logvec[num] <- logel
        num <- num + 1
    }
    tval <- 2 * (logel00 - logel)
    constraint <- as.vector(muvec %*% indicmat %*% t(nuvec))
    list(xd1 = xd1, yd1 = yd1, temp3 = temp3, mean = mean, funNPMLE = funNPMLE, 
        logel00 = logel00, logel = logel, `-2LLR` = tval, Pval = 1 - 
            pchisq(tval, df = 1), logvec = logvec, sum_muvec = sum(muvec), 
        sum_nuvec = sum(nuvec), constraint = constraint)
}