maxLikelihood <- function(data, model = c("pb", "diri"), init = NULL, maxit = 500, hess = F)
{
	p = dim(data)[2]
	
	if(model == "pb")
	{
		lhood <- function(vec)
		{
#print(vec)
			out <- pbDensity(p = p, pts = data, 
						alpha = vec[1],
						beta = vec[2:(choose(p,2)+1)])
			#print(log(out))
			#print(sum(log(out)))			
			return(-sum(log(out)))
		}
		if(is.null(init))
		{
			init <- c(1, rep(1, choose(p, 2)))
		}
	}
	else if (model == "diri")
	{
		lhood <- function(vec)
		{
#print(vec)
			out <- diriDensity(p = p, pts = data, parms = vec,
								 log = T)
#			out <- diriSpecMsr(p = p, pts = data, parms = vec)/p
#print(out)
			#print(sum(log(out)))
#			return(-sum(log(out)))
			return(-sum(out))
		}
		if(is.null(init))
		{
			init <- c(rep(1, p))
		}
	}

	#print(init)
	min <- rep(1e-5, choose(p,2)+1)
	max <- rep(1e6, length(init))
	#print(min)
	#print(lhood(init))
	#print(lhood(min))
	opt <- optim(init, lhood, method = "L-BFGS-B",lower = min, upper = max, control = list(maxit = maxit), hessian = hess)

	
	rtn <- list()
	rtn$par <- opt$par
	rtn$value <- opt$value
	rtn$convergence <- opt$convergence
	rtn$w <- data
	if(hess == F)
	{
		rtn$hessian <- 1
	}
	else
	{
		rtn$hessian <- opt$hessian
	}
	
	#likelihood contributions
	par <- opt$par
	likeContrib <- pbDensity(p = p, pts = data, alpha = par[1],
						beta = par[2:(choose(p,2)+1)])
	rtn$likeContrib <- likeContrib
	
	return(rtn)
}
