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))


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

enter image description here

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?

Thanks in advance!


1 Answer 1


I think that´s because the histograms show the marginal distributions, thus the marginal mode for each component of the matrix. The theoretical mode is the joint mode, considering all the components together. In other words, the highest probability is not achieved by optmizing individually the marginal probabilities. Sadly, I don't know a way to graphically represent the distribution's mode. But, using optimization techiniques you should be able to arrive at the correct optimum.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.