R/tbscov.R

In musto101/wilcox_R: What the package does (short line)

tbscov <-
function(x,eps=1e-3,maxiter=20,r=.45,alpha=.05){
#        Rocke's contrained s-estimator
#   returns covariance matrix only. For both locatiion and scatter, use tbs
#
#      r=.45 is the breakdown point
#      alpha=.05 is the asymptotic rejection probability.
#
if(!is.matrix(x))stop("x should be a matrix with two or more columns")
x<-elimna(x)
library(MASS)
temp<-cov.mve(x)
t1<-temp$center
s<-temp$cov
    n <- nrow(x)
    p <- ncol(x)
if(p==1)stop("x should be a matrix with two or more columns")
c1M<-cgen.bt(n,p,r,alpha,asymp=FALSE)
c1<-c1M$c1
if(c1==0)c1<-.001 #Otherwise get division by zero
M<-c1M$M
    b0 <- erho.bt(p,c1,M)
    crit <- 100
    iter <- 1
    w1d <- rep(1,n)
    w2d <- w1d
    while ((crit > eps)&(iter <= maxiter))
    {
        t.old <- t1
        s.old <- s
        wt.old <- w1d
        v.old <- w2d
        d2 <- mahalanobis(x,center=t1,cov=s)
        d <- sqrt(d2)
        k <- ksolve.bt(d,p,c1,M,b0)
        d <- d/k
        w1d <- wt.bt(d,c1,M)
        w2d <- v.bt(d,c1,M)
        t1 <- (w1d %*% x)/sum(w1d)
        s <- s*0
        for (i in 1:n)
        {
            xc <- as.vector(x[i,]-t1)
            s <- s + as.numeric(w1d[i])*(xc %o% xc)
        }
        s <- p*s/sum(w2d)
        mnorm <- sqrt(as.vector(t.old) %*% as.vector(t.old))
        snorm <- eigen(s.old)$values[1]
        crit1 <- max(abs(t1 - t.old))
#        crit <- max(crit1,crit2)
        crit <- max(abs(w1d-wt.old))/max(w1d)
        iter <- iter+1
    }
#    mnorm <- sqrt(as.vector(t1) %*% as.vector(t1))
#    snorm <- eigen(s)$values[1]
#    return(list(t1=t1,s=s))
s
}