#S-PLUS FUNCTION TO GENERATE INVERSE WISHARTS
# ----------------------------------------------------
#
# Delivers a draw from the IW distribution in k.k dimensions, with degrees of 
# freedom df and scale matrix S 
#
# For example, in the normal random sampling model with Jeffreys' prior,
#  the IW posterior for the variance matrix is sampled by the call 
#
#     sigma _  riwish(k,n-1,sighat*(n-k)) 
#
#  where sighat is the sample variance matrix based on n observations
# 
# Based on theory in P.L.Odell & A.H. Feiveson(JASA 1966 p.199-203)


riwish <- function(k,df,S)
{
   R <- diag(sqrt(2*rgamma(k,(df+k-1:k)/2)))

   R <- diag(sqrt(2*rgamma(k,(df+1-1:k)/2)))

   
   R[outer(1:k,1:k,"<")] <- rnorm (k*(k-1)/2)
   R <- t(solve(R))%*% chol(S)
   return(t(R)%*%R)
}