mixDiriDensity <- function(p = 4, k = 3, 
	pi = c(1/3, 1/3, 1/3), nu = c(2, 3, 4), 
	mu = matrix( c(1/4,1/4,1/4,1/4, 0.1,0.2,0.3,0.4, 
			0.4, 0.3, 0.2, 0.1), byrow = T,ncol = 4),
	pts = c(1/4, 1/4, 1/4, 1/4))
{
	#evaluates the density (mult by p to get spectral measure)
	#of the mixture of Dirichlet densities proposed by Boldi
	#and Davison
	
	#p is dimension of the random vector
	#k is number of components in the mixture
	#pi is the vector of weights
	#nu is the concentration parameter
	#mu is matrix of params, sum_i mu_ij = 1
	#pts are points on the simplex that density gets evaluated at
	
	if(!is.matrix(pts))
	{
		pts <- t(pts)
	}
	
	#check conditions
	checkMu <- apply(mu, 1, sum) - 1
	if(any(abs(checkMu) >  1e-6))
	{stop("mu parameters must sum to 1 for each mixture component")}
	checkCOM <- pi %*% mu - 1/p
	if(any(abs(checkCOM) >  1e-6))
	{stop("does not meet center-of-mass conditions")}
	if(abs(sum(pi) -1) > 1e-8)
	{stop("mixture weights do not sum to 1")}
	
	#evaluates the Dirichlet density for each component
	runningSum <- numeric(length = dim(pts)[1])
	for(i in seq(1,k))
	{	
		parms <- nu[i] * mu[i,]
		out <- pi[i] * diriStdDensity(p = p, parms = parms, 
				pts = pts)
		runningSum <- runningSum + out
	}
	
	return(runningSum)
}