R/discr_coords.R

In tclust: Robust Trimmed Clustering

discr_coords <- function (x, equal.weights)
{
    if (is.null (x$par$x))
        stop("Dataset not included in tclust2 object - cannot calculate discriminant coordinates.")
    
    if (missing(equal.weights))
        equal.weights <- x$par$equal.weights
    
    x.no <- x$par$x[drop = FALSE, x$cluster != 0,]
    
    p <- ncol (x.no)
    n <- nrow (x.no)
    
    W <- matrix (0, ncol = p, nrow = p)
    cl.num = x$k # sum (x$usek)
    
    if (is.null (x$cov))		# tk-means items don't have cov info
    {
        if (is.null (x$par$x))
            stop ("Cannot access the data matrix. Set store.x = TRUE.")
        x$cov <- array (dim = c (p, p, cl.num))
        for(i in 1:cl.num)
    	   x$cov[,,i] <- cov (x$par$x[x$cluster == i,])
    }
    
    for (i in 1:cl.num) {
        if (equal.weights)
            W <- W + (n - 1) / cl.num * .getsubmatrix (x$cov, i)
        else
            W <- W + (x$size[i] - 1) * .getsubmatrix (x$cov, i)
    }
    
    wm <- eigen(W, symmetric = TRUE)
    wm$values[wm$values < 1e-06] <- 1e-06
    Tm <- t (wm$vectors %*% diag (sqrt (wm$values)))
    
    #Tm <- tdecomp.hh(W)            ##  
    Tinv <- solve (Tm)
    
    S <- (n - 1) * cov (x.no)
    B <- S - W
    Z <- t (Tinv) %*% B %*% Tinv
    dc <- eigen (Z, symmetric = TRUE)
    units <- Tinv %*% dc$vectors * sqrt (n - cl.num)
    proj <- x$par$x %*% units
    
    return (proj)
}