This function samples normal inverse-wishart matrices.
Details
Consider \((Y_i, \Sigma_i) \sim MIW(M, U, \Psi, \nu)\).
Generate upper triangular factor of \(\Sigma_i = C_i C_i^T\) in the upper triangular Bartlett decomposition.
Standard normal generation: n x k matrix \(Z_i = [z_{ij} \sim N(0, 1)]\) in row-wise direction.
Lower triangular Cholesky decomposition: \(U = P P^T\)
\(A_i = M + P Z_i C_i^T\)