This function samples one matrix IW matrix.
Details
Consider \(\Sigma \sim IW(\Psi, \nu)\).
Upper triangular Bartlett decomposition: k x k matrix \(Q = [q_{ij}]\) upper triangular with
\(q_{ii}^2 \chi_{\nu - i + 1}^2\)
\(q_{ij} \sim N(0, 1)\) with i < j (upper triangular)
Lower triangular Cholesky decomposition: \(\Psi = L L^T\)
\(A = L (Q^{-1})^T\)
\(\Sigma = A A^T \sim IW(\Psi, \nu)\)