# Mode of inverse-Wishart distribution (sampled vs. calculated)

I recently investigated the sampling behavior of covariance matrices through simulation. I noticed that the mode of simulated inverse-Wishart distributed matrices somehow differs from the "theoretical" mode that can be calculated directly from the parameters of the distribution.

Suppose, I am interested in a matrix $\Sigma$ of size $p \times p$ with known entries. I denote the inverse-Wishart distribution $W^{-1}(v,S)$, where $v$ are the degrees of freedom, and $S$ is the scale matrix of size $p \times p$. Now, I understand that the mode of the inverse-Wishart is

$$\frac{S}{v+p+1} \; .$$

Therefore, simulating from $W^{-1}(v,(v+p+1)\Sigma)$ should result in a distribution of matrices with a mode that matches the entries in $\Sigma$. To my surprise, it did not!

Printed below is the code that I used to generate the matrices in R.

sim.invWishart <- function(N, v, Scale){

invSc <- solve(Scale)
p <- ncol(Scale)

invS <- rWishart(N,v,invSc)
S <- apply(invS, 3, solve)
S <- array(S,dim=c(p,p,N))

return(S)
}

df <- 10
Sigma <- matrix(c(.8,.3,.1,.3,.8,.1,.1,.1,.4),3,3)
Scale <- Sigma*(df+ncol(Sigma)+1)

N <- 50000
S <- sim.invWishart(N, df, Scale)

# mode of entries
apply(S, 1:2, function(x){d<-density(x); d$x[which.max(d$y)]})


I expected that the mode of the entries would replicate Sigma. However, the last call (mode via kernel denisities) results in the following matrix, containing the mode of each entry. I also printed a histogram to better show what I mean.

          [,1]      [,2]      [,3]
[1,] 1.1364036 0.3459420 0.1369788
[2,] 0.3459420 1.1941119 0.1168395
[3,] 0.1369788 0.1168395 0.5587053


The off-diagonal entries actually met my expectations quite well. However, the mode of the diagonal entries is always a bit off. Suprisingly, when simulating from $W^{-1}(v,v\Sigma)$ instead, the modes obtained from the simulation exactly match $\Sigma$) even though the formula for the mode suggests otherwise.

Why is that? Is my procedure flawed (e.g., simulation, calculation of the modes), or is this common behavior of the inverse-Wishart? Why does multiplying with $v+p+1$ give the wrong results, and why are they (seemingly) correct with just $v$ in that place?