A simulation-based approach to estimating the probability of an arbitrary region under a multivariate normal distribution is developed. In specific, the probability is expressed as the ratio of the unrestricted and the restricted multivariate normal density functions, where the restriction is given by the region whose probability is of interest. The density function of the restricted distribution is then estimated by using a sample generated from the Gibbs sampling algorithm.