R/concsim.R

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

concsim <-
function(n = 100, p = 2, steps = 5, gam = 0.4, runs = 20)
{
# This Splus function is used to determine when the DD
# plot separates outliers from non-outliers for various starts.
	A <- sqrt(diag(1:p))
	mbact <- 0
	fmcdct <- 0
	mbct <- 0
	madct <- 0
	dgkct <- 0
	for(i in 1:runs) {
		x <- matrix(rnorm(n * p), ncol = p, nrow = n)	
	## outliers have mean (10, 10 sqrt(2), ..., 10 sqrt(p))^T
		val <- floor(gam * n)
		tem <- 10 + 0 * 1:p
		x[1:val,  ] <- x[1:val,  ] + tem
		x <- x %*% A	#MBA
		out <- covmba(x, csteps = steps)
		center <- out$center
		cov <- out$cov
		rd2 <- mahalanobis(x, center, cov)
		if(min(rd2[1:val]) > max(rd2[(val + 1):n]))
mbact <- mbact + 1	
	#DGK
		covs <- var(x)
		mns <- apply(x, 2, mean)	## concentrate
		for(i in 1:steps) {
			md2 <- mahalanobis(x, mns, covs)
			medd2 <- median(md2)
			mns <- apply(x[md2 <= medd2,  ], 2, mean)
			covs <- var(x[md2 <= medd2,  ])
		}
		rd2 <- mahalanobis(x, mns, covs)
		if(min(rd2[1:val]) > max(rd2[(val + 1):n])) dgkct <- dgkct + 1	
	#Median Ball start
		covv <- diag(p)
		med <- apply(x, 2, median)
		md2 <- mahalanobis(x, center = med, covv)
		medd2 <- median(md2)	## get the start
		mns <- apply(x[md2 <= medd2,  ], 2, mean)
		covs <- var(x[md2 <= medd2,  ])	## concentrate
		for(i in 1:steps) {
			md2 <- mahalanobis(x, mns, covs)
			medd2 <- median(md2)
			mns <- apply(x[md2 <= medd2,  ], 2, mean)
			covs <- var(x[md2 <= medd2,  ])
		}
		rd2 <- mahalanobis(x, mns, covs)
		if(min(rd2[1:val]) > max(rd2[(val + 1):n])) mbct <- mbct + 1	
	#MAD start
		tem <- apply(x, 2, mad)^2
		covv <- diag(tem)
		md2 <- mahalanobis(x, center = med, covv)
		medd2 <- median(md2)	## get the start
		mns <- apply(x[md2 <= medd2,  ], 2, mean)
		covs <- var(x[md2 <= medd2,  ])	## concentrate
		for(i in 1:steps) {
			md2 <- mahalanobis(x, mns, covs)
			medd2 <- median(md2)
			mns <- apply(x[md2 <= medd2,  ], 2, mean)
			covs <- var(x[md2 <= medd2,  ])
		}
		rd2 <- mahalanobis(x, mns, covs)
		if(min(rd2[1:val]) > max(rd2[(val + 1):n])) madct <- madct + 1	
	#FMCD
		out <- cov.mcd(x)
		center <- out$center
		cov <- out$cov
		rd2 <- mahalanobis(x, center, cov)
		if(min(rd2[1:val]) > max(rd2[(val + 1):n]))
			fmcdct <- fmcdct + 1
	}
	list(mbact = mbact, fmcdct = fmcdct, dgkct = dgkct, mbct = mbct, madct
		 = madct)
}