"mesthub"<-
function(x, delta = 0.1)
{
# Huber's M-estimator for covariance
# x is the data matrix, delta a tuning parameter, which is usually = 0.1.
# To check for convergence the Fisher-transform of the correlation coefficient
# is used
	np <- dim(x)
	n <- np[1]
	p <- np[2]
	k <- sqrt(qchisq(1 - delta, p))
	const <- pchisq(k^2, p + 2) + ((k^2 * delta)/p)
	covstar <- var(x)
	locstar <- apply(x, 2, mean)
	fishstar <- 10^2
	step <- 1
	fish <- 10^2
	while((step <= 100) && (fish > 10^(-5))) {
		covm <- covstar
		loc <- locstar
		fishtran <- fishstar
		rd <- sqrt(mahalanobis(x, center = loc, cov = covm))
		w1 <- (rd <= k) + (k * (rd > k))/rd
		w2 <- w1^2/const
		locstar <- apply(w1 * x, 2, sum)/sum(w1)
		y <- sqrt(w2) * scale(x, center = locstar, scale = F)
		covstar <- ((n - 1) * var(y))/n
		D <- diag((diag(covstar))^(-1/2))
		corr <- D %*% covstar %*% D
		corrfish <- corr[lower.tri(corr)]
		corrfish <- (1 + corrfish)/(1 - corrfish)
		fishstar <- log(corrfish)/2
		fish <- mean(abs(fishtran - fishstar))
		step <- step + 1
	}
	if(step > 100) {
		print("Convergence failed!")
	}
	list(loc = locstar, cov = covstar, step = step)
}