#
#  Generates random matrices with a Wishart distribution
#
#  Uses the algorithm on pp. 231-232 of W J Kennedy, Jr,
#  and J E Gentle, Statistical Computing, Dekker 1980
#  taken from P L Odell and A H Feiveson, A numerical procedure
#  to generate a sample covariance matrix, J. Amer. Stat. Assoc.,
#  61 (No. 313) (1966), 199-203.
#
#  On input, n is the sample size and v is a p x p matrix-covariance 
#  matrix.
#  On output, the value returned by the function is a Wishart matrix
#  with n - p degrees of freedom and so mean (n - p)v.
#
wish <- function(n,v)
{
  p <- nrow(v)
  #  Step 1. Determine a such that v = a %*% Transpose(a)
  #  (Cholesky decomposition)
  a <- t(chol(v))
  #  Step 2.  Generate z[i][j] ~ N(0,1)
  # for i = 1, 2, ..., j and j = 2, 3, ..., p
  z <- matrix(0,p,p)
  for (j in 2:p)
    for (i in 1:j)
      z[i,j] <- rnorm(1)
  #  Step 3.  Generate y[i] ~ chi^2_{n-i}
  #  for i = 1, 2, ..., p
  y <- rchisq(p,n-i)
  #  Step 4.  Find b as in W J Kennedy, Jr, & J E Gentle, p. 231
  b <- matrix(0,p,p)
  b[1,1] <- y[1]
  for (j in 2:p) {
    b[j,j] <- y[j]
    for (i in 1:(j-1))
      b[j,j] <- b[j,j] + z[i,j]*z[i,j]
    b[1,j] <- z[1,j]*sqrt(y[1])
    }
  for (i in 2:p) {
    j <- i + 1
    while (j <= p) {
      b[i,j] <- z[i,j]*sqrt(y[i])
      k <- 1
        while (k < i) {
          b[i,j] <- b[i,j] + z[k,i]*z[k,j]
          k <- k + 1
          }
      j <- j + 1
      }
    }
  # Fill in the values of b[i,j] for i > j by symmetry
  b <- b + t(b) - diag(diag(b))
  #  Step 5.  Find s = a %*% b %*% Transpose(a)
  s <- a %*% b %*% t(a)
  return(s)
}