function(x, m)
{
	n <- nrow(x)
	k <- ncol(x)
	nk <- n * k
	av <- c(1:k)
	for(i in 1:k) {
		av[i] <- mean(x[, i])
	}
	avg <- mean(av)
	fst <- (n * var(av))/mean(diag(var(x
		)))
	fsts <- c(1:m)
	for(i in 1:m) {
		xs <- jay(n, k)
		avs <- c(1:k)
		for(j in 1:k) {
			xs[, j] <- sample(x[,
				j], n, 
				replace = T)
			avs[j] <- mean(xs[, 
				j])
		}
		avsg <- mean(avs)
		devs <- c(1:k)
		for(j in 1:k) {
			devs[j] <- n * (avs[
				j] - av[j])^
				2
		}
		fsts[i] <- (sum(devs) - nk * (
			avg - avsg)^2)/((k - 
			1) * mean(diag(var(
			xs))))
	}
	prop <- 0 * jay(m, 1)
	for(i in 1:m) {
		if(fsts[i] >= fst)
			prop[i] <- 1
	}
	pvalue <- mean(prop)
	res <- list(fst, pvalue)
	names(res) <- c(" F for ANOVA", 
		"P-Value")
	res
}
