HPD.gamma=function(n,ybar,guesses,C,niter=20){
# 
#
#  This function finds values new1 and new2 such that a
#  Gamma(n,n*ybar) density equals C.  It then computes the probability, p1,
#  to the left of new1 and the probability, p2, to the right of new2.
#  
#  Guesses is a vector of two elements containing initial guesses for
#  new1 and new2.
#
#  The quantity niter is the number of iterations used
#  by Newton's method to approximate new1 and new2.
#  
#
alpha=n
beta=n*ybar
con=alpha*log(beta)-lgamma(alpha)-log(C)
new1=guesses[1]
new2=guesses[2]
for(i in 1:niter){
new1=new1-(con+(alpha-1)*log(new1)-beta*new1)/((alpha-1)/new1-beta)
}
for(i in 1:niter){
new2=new2-(con+(alpha-1)*log(new2)-beta*new2)/((alpha-1)/new2-beta)
}
p1=pgamma(new1,alpha,beta)
p2=1-pgamma(new2,alpha,beta)
c(new1,new2,p1,p2,p1+p2)
}