#
#  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.
#  On input, p is the dimension of the matrix, n the sample size
#  On output, w is the Wishart matrix and invw its inverse
#
wish <- function(p,n,v)
{
  for (i in 1:p)
    for (j in 1:p)
      if (i < j) v[i,j] <- v[j,i]
  #  Step 1. Determine a such that v = a*+Transpose(a)
  #  d consists of eigenvalues, evec of eigenvectors
  a <- matrix(0,p,p)
  for (i in 1:p)
    for (j in 1:p)
      a[i,j] <- eigen(v)$vectors[i,j]*sqrt(eigen(v)$values[j])
  #  Step 2.  Generate z[i][j] ~ N(0,1)
  #  for i <= j = 1, 2, ..., p
  z <- matrix(0,p,p)
  for (j in 1: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[i])
    for (i in 1:(j-1)) {
      b[i,j] <- z[i,j]*sqrt(y[i])
      for (k in 1:(i-1))
        b[i,j] <- b[i,j] + z[k,i]*z[k,j]
      }
    for (i in 1:j)
      b[i,j] <- b[j,i]
  }
  #  Step 5.  Find s = a*+b*+Transpose(a)
  #  First intermed = b*+Transpose(a)
  intermed <- matrix(0,p,p)
  for (i in 1:p)
    for (j in 1:p)
      for (k in 1:p)
        intermed[i,j] <- intermed[i,j] + b[i,k]*a[j,k]
  #  Now s = a*+intermed
  s <- matrix(0,p,p)
  for (i in 1:p)
  for (j in 1:p) {
    s[i,j] <- 0
    for (k in 1:p)
      s[i,j] <- s[i,j] + a[i,k]*intermed[k,j]
    }
  #   Now ensure symmetry
  for (i in 1:p)
    for (j in 1:(i-1)) {
      temp <- (s[i,j]+s[j,i])/2
      s[i,j] <- temp
      s[j,i] <- temp
    }
  #  Now s is a Wishart matrix (aside
  #  from the factor 1/n) with n-p d.f.
  #  (W J Kennedy, Jr, & J E Gentle, p. 232)
  #
  #  Step 6.  Find inverse invs of s
  #  First determine a such that
  #  s = a*+Transpose(a)
  #  Resymmetrize s
  for (i in 1:p)
    for (j in 1:(i-1))
      s[j,i] <- s[i,j]
  w <- s
  return(b)
}