R/corrsim.R

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

corrsim <-
function(n = 100, p = 3, eps = 0.4, nruns = 100, type = 1)
{
#For R, first type "library(lqs)" before using this function
# This function generates 100 n by p matrices x.
# The output is the 100 sample correlations between the MDi and RDi
# RDi uses covmba for type = 1, rmba for type = 2, cov.mcd for type = 3
# mahalanobis gives squared Maha distances
	corrs <- 1:nruns
	for(i in 1:nruns) {
		wt <- 0 * (1:n)
		x <- matrix(rnorm(n * p), ncol = p, nrow = n)	
	#The following 3 commands make x elliptically contoured.
#zu <- runif(n)
#x[zu < eps,] <- x[zu < eps,]*5
#x <- x^2
# To make marginals of x lognormal, use
#x <- exp(x)
		center <- apply(x, 2, mean)
		cov <- var(x)
		md2 <- mahalanobis(x, center, cov)
		if(type == 1) {
			out <- covmba(x)
		}
		if(type == 2) {
			out <- rmba(x)
		}
		if(type == 3) {
			out <- cov.mcd(x)
		}
		center <- out$center
		cov <- out$cov
		rd2 <- mahalanobis(x, center, cov)	
	# need square roots for the usual distances
		md <- sqrt(md2)
		rd <- sqrt(rd2)
		const <- sqrt(qchisq(0.5, p))/median(rd)
		rd <- const * rd	
	# wt[rd < sqrt(qchisq(0.975, p))] <- 1
#  corrs[i] <- cor(md[wt > 0], rd[wt > 0])}
		corrs[i] <- cor(md, rd)
	}
	cmean <- mean(corrs)
	cmin <- min(corrs)
	clt95 <- sum(corrs < 0.95)
	clt80 <- sum(corrs < 0.8)
	list(cmean = cmean, cmin = cmin, clt95 = clt95, clt80 = clt80,
		corrs = corrs)
}