require(MASS) # Needed to generate multivariate random normals
require(stats) # For factor analysis; but if this isn't loaded
			   # you have bigger problems!

factor.from.random.corr = function(d, n=1000, maxcor=2/(d-1), 
								   mincor=0, output="variance") {
	# Get a random correlation matrix
	sigma = make.corr.matrix(d,maxcor=maxcor,mincor=mincor,
	        				 diag.dom=FALSE)
	# Multiply by the standard variance of 15^2 to go from the
	# correlation matrix to the covariance matrix
	varcovar = sigma*15*15
	# set the means to 100
	mu = rep(100,d)
	# Make up n random Gaussian vectors
	x = mvrnorm(n,mu,varcovar)
	# Do maximum-likelihood factor analysis with 1 factor
	my.fa = factanal(x,1)
	v = as.vector(my.fa$loadings)
	switch(output,
		# Extract the fraction of variance described by the factor
		variance = sum(v^2)/d,
		loadings = v,
		everything = list(loadings=v,sigma=sigma,varexp=sum(v^2)/d,data=x),
		"I can't tell what you want me to output")
}


make.corr.matrix = function(d, maxcor=2/(d-1), mincor=0, 
							diag.dom=FALSE) {
	# Create a random symmetric positive-definite matrix with
	# all non-negative entries to serve as a correlation matrix
	# It's easy to get something symmetric with non-negative
	# entries...
	# Use a rejection procedure: keep drawing from the easy
	# distribution until you get something positive-definite
	if (diag.dom) {
		# To make a diagonally-dominated matrix, create a test
		testf = function(m) {
			# A symmetric matrix which is diagonally dominated
			# is positive symmetric.  "Diagonally dominated"
			# means (1) the diagonal entries are positive, and
			# (2) for each row, the sum of the absolute values of
			# the non-diagonal entries is less than the diagonal
			# entry.  Since this is a correlation matrix, all of
			# the diagonal entries will be 1, taking care of the
			# first criterion, and simplifying the second.
			diag(m) = 0 # Doesn't change the outside matrix and
						# simplifies doing the row-sums
			diag.dominated = TRUE
			# Test each rom for diagonal domination
			for  (i in 1:d) {
				# Loops are slow in r, but d is small here...
				rowsum = sum(m[i,])
				# If any row fails, the whole matrix fails
				if (rowsum >= 1) {
					diag.dominated = FALSE
				}
			}
			return(diag.dominated)
		} }
	else {
		# Create any matrix so long as it's positive-definite
		# so define a test just for that
		testf = function(m) {
			# chol(m) attempts to perform the Cholesky
			# decomposition of m, which only makes sense if m
			# is symmetric and positive-definite.  It gives
			# an error if m is not positive definite (but doesn't
			# check symmetry).  try(expr,silent=TRUE) returns
			# the value of the expression, or nothing if it
			# produced an error.  So try(chol(m),silent=TRUE)
			# either gives us a matrix, or nothing.  So
			# wrapping the whole thing in is.matrix() tells
			# us whether m was positive definite or not.
			# Admittedly, it does waste some time by actually
			# doing the Cholesky decomp. on m if it can...
			p.d = is.matrix(try(chol(m),silent=TRUE))
			return(p.d)
		} }
	pos.def = FALSE
	while (!pos.def) {
			sigma = make.symm.matrix(d, maxentry=maxcor, 
									 minentry=mincor)
			diag(sigma) = 1
			pos.def = testf(sigma)
	}
	return(sigma)
}

# Create a symmetric matrix with uniformly-distributed positive
# entries, and zeroes on the diagonal.
make.symm.matrix = function(d,minentry=0,maxentry=1) {
	# Number of distinct components = no. in upper triangular
	# part above diagonal
	# (d-1) + (d-2) + ... + 1= (d-1)d/2
	# Start with an all-zero matrix
	m = matrix(data=0,nrow=d,ncol=d)
	# Create the entries
	entries = runif(d*(d-1)/2,min=minentry,max=maxentry)
	for (i in 1:(d-1)) {
		for (j in (i+1):d) {
			m[i,j] = entries[1]
			m[j,i] = entries[1]
			entries = entries[-1] # Remove the first element from
								  # the list
		}
	}
	return(m)
}