library(MCMCpack)

# Define a function that computes the moments of an inverse Wishart distribution
# given a scale matrix V and degrees of freedom v
inv.wishart.moments <- function(V, v) {
 # some sanity checks upfront
 stopifnot(is.matrix(V), nrow(V) == ncol(V))
 p <- ncol(V)
 stopifnot(v > p - 1)

 # compute the mean matrix
 if (v-p-1 != 0) {
  mean.matrix <- V / (v - p - 1)
 } else {
  mean.matrix <- NA
 }

 # compute the variances
 if ( (v-p-3) != 0 && (v-p-1) != 0 && (v-p) !=0 ) {
  var.matrix <- matrix(NA, ncol=p, nrow=p)
  for (i in 1:p) {
   for (j in 1:p) {
    var.matrix[i,j] <- ( (v-p+1)*V[i,j]^2+(v-p-1)*V[i,i]*V[j,j] ) / ((v-p)*(v-p-1)^2*(v-p-3))
   }
  }
  inv.var.matrix <- solve(var.matrix)
 } else {
  var.matrix <- NA
  inv.var.matrix <- NA
 }

 # return the results
 return(list(mean.matrix=mean.matrix, var.matrix=var.matrix, inv.var.matrix=inv.var.matrix))
}

# Define a function that simulates draws from a Wishart and inverse Wishart distribution
# Draw 10000 samples and compute the means and variances.
empirical.wishart.moments <- function(V, v) {
 # some sanity checks
 stopifnot(is.matrix(V), nrow(V) == ncol(V))
 p <- ncol(V)

 a <- array(dim=c(10000,2,2))

 # sample the inverse Wishart first
 for(i in 1:10000) { a[i,,] <- riwish(v, V) }
 m11 <- mean(a[,1,1]); 
 m12 <- mean(a[,1,2]);
 m22 <- mean(a[,2,2]);
 m21 <- mean(a[,2,1]);
 v11 <- var(a[,1,1]); 
 v12 <- var(a[,1,2]);
 v22 <- var(a[,2,2]);
 v21 <- var(a[,2,1]);
 m1 <- matrix(c(m11, m12, m21, m22), nrow=2, byrow=TRUE)
 v1 <- matrix(c(v11, v12, v21, v22), nrow=2, byrow=TRUE)

 # sample the corresponding Wishart,
 # Notice the parameterization of m=v-p-1
 for(i in 1:10000) { a[i,,] <- rwish(v - p - 1, solve(V)) }
 m11 <- mean(a[,1,1]); 
 m12 <- mean(a[,1,2]);
 m22 <- mean(a[,2,2]);
 m21 <- mean(a[,2,1]);
 v11 <- var(a[,1,1]); 
 v12 <- var(a[,1,2]);
 v22 <- var(a[,2,2]);
 v21 <- var(a[,2,1]);
 m2 <- matrix(c(m11, m12, m21, m22), nrow=2, byrow=TRUE)
 v2 <- matrix(c(v11, v12, v21, v22), nrow=2, byrow=TRUE)

 return(list(inv.wishart.means=m1, inv.wishart.variances=v1, wishart.means=solve(m2), wishart.variances=solve(v2)))

}

# create a sequence of different Wishart distributions with the same mean
# but different degrees of freedom. Note that the variances decrease with
# increasing df
# Thus, if we want a large variance (i.e. as an uninformative prior), we must
# choose small degrees of freedom.
for(v in c(5,6,10,50,100)) {
  V <- matrix(c(1*(v-3), 0.5*(v-3), 0.5*(v-3), 1*(v-3)), nrow=2, byrow=TRUE)
  print(v, digits=2)
  print(inv.wishart.moments(V, v)[c("mean.matrix", "var.matrix")], digits=2)
  print(empirical.wishart.moments(V, v)[c("inv.wishart.means", "inv.wishart.variances")], digits=2)
}

