condDensApprox <- function(p, z, model = "pairwiseBeta", par, 
					predLoc = p, numInitPts = 201)
{
	
	# z is the entire vector of observations (observed & missing)
	# p is the dimension (including missing) 
	# model is the angular density model
	# par is the parameter vector needed for the ang dens model
	# predLoc is the location assumed missing
	# numInitPts is number of points chosen on the simplex
	
	#zNot represents the unobserved value
	#w's are the points on the spectral density
	zObs <- z[-predLoc]
	zMatrix <- matrix(rep(z, numInitPts), nrow = numInitPts, 
					byrow = T)
	zObsMtx <- zMatrix
	zObsMtx[, predLoc] <- NA	

	#creates the set of points on the simplex first (the line 
	#is defined by the (p-1) known components
	wNot <- seq(0, 1, length = numInitPts)
	divisor <- matrix(rep(sum(zObs)/(1-wNot), p), ncol = p)
	wSupp <- zObsMtx/divisor
	wSupp[, predLoc] <- wNot
	w <- wSupp
	
	#gives corresponding radii
	if(predLoc != 1)
	rad <- z[1]/w[,1]
	else
	rad <- z[2]/w[,2]

	#gives corresponding values of zNot
	zNot <- wNot*rad

	#gets corresponding values of xNot
	xNot <- gInvFun(zNot)
	
	#fills in some lower values of xNot
	xNotLow <- seq(xNot[1], xNot[2], 2)
	xNot <- c(xNotLow, xNot[-1])
	
	#backtransforms to get zNot values
	zNotLow <- gFun(xNotLow)
	zNot <- c(zNot[1], zNotLow[-1], zNot[-1])
	
	#gets an odd number of nodes
	if(length(xNot) %% 2 == 0)
	{
		zNot <- zNot[-length(xNotLow)]
		xNot <- xNot[-length(xNotLow)]
	}
	numPts <- length(zNot)
	
	#final z's and w's to be used
	zMatrix <- matrix(rep(z, numPts), nrow = numPts, 
					byrow = T)
	zObsMtx <- zMatrix
	zObsMtx[, predLoc] <- zNot
	rad <- apply(zObsMtx, 1, sum)
	w <- zObsMtx/rad

	
	#gets values of angular measure at these w locations
	#used only for returning, integration done below
    #gives the values of the angular measure at the above points
	if(model == "pairwiseBeta")
	{
		alpha <- par[1]
		beta <- par[-1]
		smVals <- pbSpecMsr(p, alpha, beta, w) 
	}
	else if(model == "mixedDirichlet")
	{
		K <- length(par)/(2+p)
		pi <- par[1:K]
		nu <- par[(K+1):(2*K)]
		mu <- matrix(par[(2*K+1):length(par)], ncol = p, nrow = K)
		#mult below by p to get spectral measure
		smVals <- p * mixDiriDensity(p, K, pi, nu, mu, pts = w)
	}
	else
	{
		stop("Invalid angular density model")
	}

	#gets (unnormalized for now) density values for returning
	unDensVal <- rad^(-p-1) * smVals
	
	
	############################################################
	# To test the accuracy using the binning procedure we need to
	# estimate the approximate conditional distribution -- that 
	# is, we need a "running integral" of the density.  We 
	# estimate the integral by using adaptive quadrature at the
	# lower tail (slow but accurate), and using a Simpson's
	# for the interior (faster)
	############################################################

	###########
    #lower tail
    ###########
    
    #function to be integrated using adaptive quad procedure in R
    #takes a vector of zNot's, gets a vector of radii, gets
    #corresponding w's, returns the unnormalized density value
	#this depends on the angular measure model
	if(model == "pairwiseBeta")
	{
		alpha <- par[1]
		beta <- par[-1]
		
	    unDensValEval <- function(zNot)
    	{
			zTemp <- matrix(rep(z, length(zNot)), byrow = T, 
						ncol = p)      
			zTemp[, predLoc] <- zNot
    		rad <- apply(zTemp, 1, sum)
    		w <- zTemp/rad
	    	smVal <- pbSpecMsr(p, alpha, beta, w)
    		unDensVal <- rad^(-p-1)*smVal
    		return(unDensVal)
    	}
    }
	else if(model == "mixedDirichlet")
	{
		K <- length(par)/(2+p)
		pi <- par[1:K]
		nu <- par[(K+1):(2*K)]
		mu <- matrix(par[(2*K+1):length(par)], ncol = p, nrow = K)
		
		unDensValEval <- function(zNot)
		{
			zTemp <- matrix(rep(z, length(zNot)), byrow = T, 
						ncol = p)      
			zTemp[, predLoc] <- zNot
    		rad <- apply(zTemp, 1, sum)
    		w <- zTemp/rad
	    	smVal <- p * mixDiriDensity(p, K, pi, nu, mu, pts = w)
    		unDensVal <- rad^(-p-1)*smVal
    		return(unDensVal)
		}
	}
    else
    {
    	stop("invalid angular density model in lower tail")
    }

	#performs the integration from zNot[1] = 0 to zNot[2] 
	#(the first non-zero value of zNot)
    lowTail <- integrate(unDensValEval, 0, zNot[2])$value
    

	###########
    #everything but lower tail
    ###########

    #expands the zNot vector to include midpoints

	#quick-and-dirty fix of the upper tail
	zNotExpand <- zNot
	zNotExpand[length(zNotExpand)] <- 2 * 
					zNotExpand[length(zNotExpand) - 1]

	#adds midpoints
    midPoints <- (zNotExpand[2:(length(zNotExpand)-1)] + 
    					zNotExpand[3:length(zNotExpand)])/2
	zNotExpand <- sort(c(zNotExpand[-1], midPoints))
    lenExpand <- length(zNotExpand)
    
	zObsMtxExpand <- matrix(rep(z, lenExpand), nrow = lenExpand, 
						byrow = T)    
	zObsMtxExpand[, predLoc] <- zNotExpand
	
	#gets angular components
	radExpand <- apply(zObsMtxExpand, 1, sum)
	wExpand <- zObsMtxExpand/radExpand
	
	#gets xNot values for later use in conditional expectation
	xNotExpand <- gInvFun(zNotExpand)

    #gives values of angular measure at the above wExpand's
	if(model == "pairwiseBeta")
	{
		alpha <- par[1]
		beta <- par[-1]
		smValsExpand <- pbSpecMsr(p, alpha, beta, wExpand) 
	}
	else if(model == "mixedDirichlet")
	{
		K <- length(par)/(2+p)
		pi <- par[1:K]
		nu <- par[(K+1):(2*K)]
		mu <- matrix(par[(2*K+1):length(par)], ncol = p, nrow = K)
		#mult below by p to get spectral measure
		smValsExpand <- p * mixDiriDensity(p, K, pi, nu, mu,
						pts = wExpand)
	}
	else
	{
		stop("Invalid angular density model")
	}

	#gives the approximate values for the unnormalized density
	#at zNot defined above
	#dens \approx 1/rad^(p+1)*smVals
	unDensValExpand <- radExpand^(-p-1) * smValsExpand
	unDensValExpand[is.nan(unDensValExpand)] <- 0

	#normalizes by using a modified simpson's rule to approximate
	#the integral--returns the cumulative integral for use
	#in finding quantiles, etc.
	
	#estimates the integral using Simpson's Rule
	outSimp <- unequalSimpson(zNotExpand, unDensValExpand)
	traps <- c(0, lowTail, outSimp)
	cumTraps <- cumsum(traps)
	normalizer <- cumTraps[length(cumTraps)]
	zDensVal <- unDensVal/normalizer
	cumTraps <- cumTraps/normalizer		


	#does the conversion to get the density values on the x-scale
	gPrime <- gDerFun(xNot)
	xDensVal <- zDensVal * gPrime
	gPrimeExpand <- gDerFun(xNotExpand)
	xDensValExpand <- unDensValExpand/normalizer * gPrimeExpand

	
	############
	# conditional expectation
	############
	#lower tail -- A LOT of work so that the lower tail 
	#can be integrated

	intFun <- function(x)
	{
		#get corresponding z Vector
		zNot <- gFun(x)
		zMtx <- matrix(rep(z, 1, each = length(x)), 
					ncol = length(z))
		zMtx[, predLoc] <- zNot
		rowSum <- apply(zMtx, 1, sum)
		wMtx <- zMtx/rowSum
		
		#get value of cond Dens on z scale
		if(model == "pairwiseBeta")
		{
			#get values of spectral measure
			smVal <- pbSpecMsr(p, par[1], par[-1], wMtx)
			#mult by radius to get f(z)
			unDensVal <- rowSum^(-p-1)*smVal
			zDensVal <- unDensVal/normalizer
		}
		else if(model == "mixedDirichlet")
		{
			K <- length(par)/(2+p)
			pi <- par[1:K]
			nu <- par[(K+1):(2*K)]
			mu <- matrix(par[(2*K+1):length(par)], ncol = p, 
						nrow = K)
			smVal <- p * mixDiriDensity(p, K, pi, nu, mu, 
						pts = wMtx)
			unDensVal <- rowSum^(-p-1)*smVal
			zDensVal <- unDensVal/normalizer
		}	

	
		#function to eval is x * f(z) * gPrime(x)
		val <- x * zDensVal * gDerFun(x)
	
		return(val)
	}
	
	lowTailCE <- integrate(intFun, xNot[1], xNot[2], 
					stop.on.error = F)$value
	
	#rest of it done by simpson's rule
	#quick fix of upper tail
	xNotInt <- xNotExpand
	xDensValInt <- xDensValExpand
	xNotInt[length(xNotInt)] <- 2 * xNotInt[length(xNotInt) - 1]
	xDensValInt[length(xDensValInt)] <- 0
	simpOutCE <- unequalSimpson(xNotInt, xNotInt * xDensValInt)
	
	xConditionalExp <- sum(c(lowTailCE, simpOutCE))	

	out <- list()
	out$zNot <- zNot
	out$zDensVal <- zDensVal
	out$cumTraps <- cumTraps 
	out$w <- w
	out$z <- w*rad
	out$smVals <- smVals
	out$zActual <- z
	out$predLoc <- predLoc
	out$xNot <- xNot
	out$xDensVal <- xDensVal
	out$xConditionalExp <- xConditionalExp
		
	return(out)
}
