#GENERATE WISHART------------------------------------------------------------
"rwish" <- function(s,nu,Cov)
{
#sxs Wishart matrix, nu degree of freedom, var/covar Cov based on
#P.L.Odell & A.H. Feiveson(JASA 1966 p.199-203). Returns w=(RU)'RU
#where Cov=U'U (U is upper triang) and where upper-tri R is
# R_ij~N(0,1), i<j ; (R_ii)^2~Chisq(nu-s+i)
if (nu<=s-1) stop ("Wishart algorithm requires nu>s-1")
R<- diag(sqrt(2*rgamma(s,(nu+1 - 1:s)/2)))
R[outer(1:s, 1:s,  "<")] <- rnorm (s*(s-1)/2)
R <- R%*% chol(Cov)
return(t(R)%*%R)
} 
#GENERATE INVERSE WISHART----------------------------------------------------
"riwish" <- function(s,df,Prec)
{
#sxs Inverse Wishart matrix, df degree of freedom, precision matrix
#Prec. Distribution of W^{-1} for Wishart W with nu=df+s-1 degree of
# freedoom, covar martix Prec^{-1}.
# NOTE  mean of riwish is proportional to Prec
if (df<=0) stop ("Inverse Wishart algorithm requires df>0")
R <- diag(sqrt(2*rgamma(s,(df + s  - 1:s)/2)))
R[outer(1:s, 1:s,  "<")] <- rnorm (s*(s-1)/2)
S <- t(solve(R))%*% chol(Prec)
return(t(S)%*%S)
}