### VERSION OF AUG/06
### cleaned up, commented, checked

igdens <- function(z, a, b) {
### returns inverse-gamma density function
tmp <- rep(0,length(z))
tmp[z>0] <- exp(-(a+1)*log(z[z>0]) - b/(z[z>0]))
tmp
} 

analyze <- function(alpha, beta, tau2, lam2, sig2,
                   sd.alpha, sd.beta, 
                   scl.tau, shp.tau,
                   scl.lam, shp.lam,
                   scl.sig, shp.sig, GRLEN=5000, CHCK=T) {

### takes as inputs value of
### theta=(alpha[1:2],beta[1:3],tau2,lam2,sig2),
### and hyperparameters sd.alpha[1:2], sd.beta[1:3],
### scl and shp for tau^2, lam^2, sig^2  
### returns conditional prior dens of phi[3] given phi[-3]=phi[1:2,4:8]
### on grid of phi[3] values, for the specified phi[-3] values
  
  ### compute parameter value in transparent parameterization
  phi <- c(alpha,
           beta[3], beta[1]+alpha[1]*beta[2], alpha[2]*beta[2]+beta[3],
           tau2+lam2,
           beta[2]*lam2,
           sig2+beta[2]^2*lam2)

  ### check coding of inverse mapping
  ### should return rep(0,8)
  if (CHCK) {
    print("check")
    print(c(alpha, beta, tau2, lam2, sig2) -
          c(phi[1:2],
           phi[4]-(phi[1]/phi[2])*(phi[5]-phi[3]),
           (phi[5]-phi[3])/phi[2],
           phi[3],
           phi[6]-phi[2]*phi[7]/(phi[5]-phi[3]),
           phi[2]*phi[7]/(phi[5]-phi[3]),
           phi[8]-(phi[7]/phi[2])*(phi[5]-phi[3]) ) ) }

  ### determine support (lft,rht) of desired conditional
  aa <- abs(phi[2]*phi[7])/phi[6]; bb <- abs(phi[2]/phi[7])*phi[8]
  if ( (phi[2]*phi[7]) > 0) {
    lft <- phi[5]-bb; rht <- phi[5]-aa }
  if ( (phi[2]*phi[7]) < 0) {
    lft <- phi[5]+aa; rht <- phi[5]+bb }

  ### is support wide relative to prior, truncate for numerical reasons
  if (rht >    5*sd.beta[3]) {rht <-    5*sd.beta[3]}
  if (lft < (-5)*sd.beta[3]) {lft <- (-5)*sd.beta[3]}

  ### set up grid for desired conditional
  gr <- seq(from=lft, to=rht,length=GRLEN)

  ### beta[3] term
  trma <- dnorm(gr, sd=sd.beta[3])
  ### beta[2] term
  trmb <- dnorm((phi[5]-gr)/phi[2], sd=sd.beta[2])
  ### beta[1] term
  trmc <- dnorm(phi[4]-(phi[1]/phi[2])*(phi[5]-gr), sd=sd.beta[1])
  ### lam2 term
  trmd <- igdens(phi[7]*phi[2]/(phi[5]-gr), shp.lam, scl.lam)
  ### tau2 term
  trme <- igdens(phi[6]-phi[7]*phi[2]/(phi[5]-gr), shp.tau, scl.tau)
  ### sig2 term
  trmf <- igdens(phi[8]-(phi[7]/phi[2])*(phi[5]-gr), shp.sig, scl.sig) 
  ### Jacobian term
  trmg <- 1/abs(phi[5]-gr)

  ### combine terms and normalize to density
  ans <- log(trma)+log(trmb)+log(trmc)+log(trmd)+
         log(trme)+log(trmf)+log(trmg)
  ans <-ans-max(ans); ans <- exp(ans)
  ans <- (ans/sum(ans))/(gr[2]-gr[1])

### return conditional density, also mean, truncation points  
list(gr=gr, pstdens=ans, 
     pstmn=sum(gr*ans)/sum(ans), lft=lft, rht=rht)
}