This paper provides a framework for estimating the mean and variance of a highdimensional normal density. The main setting considered is a fixed number of vector following a high-dimensional normal distribution with unknown mean and diagonal covariance matrix. The diagonal covariance matrix can be known or unknown. If the covariance matrix is unknown, the sample size can be as small as 2. The proposed estimator is based on the idea that the unobserved pairs of mean and variance for each dimension are drawn from an unknown bivariate distribution, which we model as a mixture of normal-inverse gammas. The mixture of normal-inverse gamma distributions provides advantages over more traditional empirical Bayes methods, which are based on a normal-normal model. When fitting a mixture model, we are essentially clustering the unobserved mean and variance pairs for each dimension into different groups, with each group having a different normal-inverse gamma distribution. The proposed estimator of each mean is the posterior mean of shrinkage estimates, each of which shrinks a sample mean towards a different component of the mixture distribution. Similarly, the proposed estimator of variance has an analogous interpretation in terms of sample variances and components of the mixture distribution. If diagonal covariance matrix is known, then the sample size can be as small as 1, and we treat the pairs of known variance and unknown mean for each dimension as random observations coming from a flexible mixture of normal-inverse gamma distributions.