#-----------------------------------------
# Atypicality index for compositional data
# Based on Aitchison(1986, 2003)
#
# January 25, 2006
# Developed by Hiroyoshi Arai
#-----------------------------------------

atypicality <- function(x, z=NULL) {

	d <- ncol(x)-1
	x <- x/apply(x, 1, sum)
	Y <- data.frame(log(x[ ,1:d]/x[ ,d+1]))

	if (!is.null(z)) {
		if (is.vector(z)) {
			z <- matrix(z, nrow=1)
		}
		z <- data.frame(z)
		if (ncol(x) != ncol(z)) {
			stop("The number of columns is different between x and z.")
		}
		z <- z/apply(z, 1, sum)
		y <- data.frame(log(z[ ,1:d]/z[ ,d+1]))
		Ymu <- mean(Y)
		Ysigma <- var(Y)
		N <- nrow(x)
	} else {
		y <- Y
		N <- nrow(x)-1
	}

	result <- data.frame(matrix(NA, nrow=nrow(y), ncol=2))
	colnames(result) <- c("q(y)", "atypicality")
	n <- N-1

	for (i in 1:nrow(y)) {
		if (is.null(z)) {
			Ymu <- mean(Y[-i, ])
			Ysigma <- var(Y[-i, ])
		}
		qy <- 1/(1+1/N)*as.matrix(y[i, ]-Ymu)%*%solve(Ysigma)%*%t(y[i, ]-Ymu)
		result[i, ] <- c(qy, pbeta(qy/(qy+n), d/2, (n-d+1)/2))
	}

	return(result)
}