pbDensity <- function(p = 4, alpha = 1, beta = c(2,2,2,2,2,2),
						pts = c(.25, .25, .25, .25))
{	
#############
#NOTE:  OCT 13, 2008
#WHILE WORKING ON PREDICITION PROBLEM, THIS FUNCTION WAS UPDATED TO BETTER
#WORK WITH MAX-LIKELIHOOD FUNCTION.  YOU CAN SEE UPDATES IN THE PREDEXTREMES
#DIRECTORY
##############


	#evaluates the pairwise beta density at a given point
	#on the simplex with given pairwise parameters for a 
	#given dimension p
	
	if(is.vector(pts))
	{
		pts <- t(pts)
	}

	count <- 0
	ptsPair <- matrix(ncol = choose(p,2), nrow = dim(pts)[1])
	ptsTop  <- matrix(ncol = choose(p,2), nrow = dim(pts)[1])
	
	for(i in seq(1,(p-1)))
	{
		for(j in seq(i+1,p))
		{
			count <- count + 1
			ptsPair[,count] <- pts[,i] + pts[,j]
			ptsTop[,count]  <- pts[,i]
		}
	}
	
	a <- 2*alpha
	b <- alpha*(p-2)-p+3
	
	one   <- (ptsPair)^(a-1)*(1-ptsPair)^(b-1)
	two   <- matrix(rep(exp(lgamma(2*beta) - 2*lgamma(beta)), 
						dim(pts)[1]), ncol = choose(p,2), 
						byrow = T)				
	three <- t((t(ptsTop/ptsPair))^(beta - 1))
	four  <- t((t(1 - ptsTop/ptsPair))^(beta - 1))
	const <- 2*gamma(p-2)/(sqrt(p) * p * (p-1)) *
			exp(lgamma(a+b+p-1) - lgamma(a+1) -
				lgamma(b+p-3) - log(a+b+p-2))

	h     <- const*one*two*three*four

	return(apply(h,1,sum))
}	
		
