Mixture of multivariate normal inverse Gaussian (MNIG) distributions could be applied to clustering data with skewness and heavy tails. In this talk, I will discuss the process that estimates the MNIG parameters under a Bayesian framework via a Gibbs scheme. I will discuss two novel approaches to simulate generalized inverse Gaussian (GIG) random numbers and matrix generalized inverse Gaussian (MGIG) random matrices, which play important role in estimating the MNIG parameters. Our algorithm will be applied to both simulated and real data.