R/cmba2.R

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

cmba2 <-
function(x, csteps = 5, ii = 1)
{
# gets the covmba estimator using 98, 95, 90, 80, 70, 60 and 50% trimming
	n <- dim(x)[1]
	p <- dim(x)[2]
	mds <- matrix(nrow = n, ncol = 8, 0)	##get the DGK estimator
	covs <- var(x)
	mns <- apply(x, 2, mean)
	cmd <- sqrt(mahalanobis(x, mns, covs))	## concentrate
	for(i in 1:csteps) {
		md2 <- mahalanobis(x, mns, covs)
		medd2 <- median(md2)
		mns <- apply(x[md2 <= medd2,  ], 2, mean)
		covs <- var(x[md2 <= medd2,  ])
	}
	mds[, 8] <- sqrt(mahalanobis(x, mns, covs))
	covb <- covs
	mnb <- mns	##get the square root of det(covb)
	critb <- prod(diag(chol(covb)))	##get the resistant estimator
	covv <- diag(p)
	med <- apply(x, 2, median)
	md2 <- mahalanobis(x, center = med, covv)
	smd2 <- sort(md2)
	val <- p + 3
	tem <- 1:7
	tem[1] <- smd2[val + floor(0.02 * n)]
	tem[2] <- smd2[val + floor(0.05 * n)]
	tem[3] <- smd2[val + floor(0.1 * n)]
	tem[4] <- smd2[val + floor(0.2 * n)]
	tem[5] <- smd2[val + floor(0.3 * n)]
	tem[6] <- smd2[val + floor(0.4 * n)]
	tem[7] <- median(md2)
	medd2 <- tem[7]
	for(j in ii:7) {
## get the start
		val2 <- tem[j]
		mns <- apply(x[md2 <= val2,  ], 2, mean)
		covs <- var(x[md2 <= val2,  ])	## concentrate
		for(i in 1:csteps) {
			md2 <- mahalanobis(x, mns, covs)
			medd2 <- median(md2)
			mns <- apply(x[md2 <= medd2,  ], 2, mean)
			covs <- var(x[md2 <= medd2,  ])
		}
		mds[, j] <- sqrt(mahalanobis(x, mns, covs))
		plot(cmd, mds[, j])
		identify(cmd, mds[, j])
		crit <- prod(diag(chol(covs)))
		if(crit < critb) {
			critb <- crit
			covb <- covs
			mnb <- mns
		}
	}
	pairs(mds)	##scale for better performance at MVN
	rd2 <- mahalanobis(x, mnb, covb)
	const <- median(rd2)/(qchisq(0.5, p))
	covb <- const * covb
	list(center = mnb, cov = covb, mds = mds)
}