The empirical Bayes $g$-modeling approach via the nonparametric maximum likelihood estimator (NPMLE) is widely used for large-scale estimation and inference in the normal means problem, yet theoretical guarantees for uncertainty quantification remain scarce. A key obstacle is that the NPMLE of the mixing distribution is necessarily discrete, which yields discrete posterior credible sets and a deconvolution rate that is logarithmic. We address both limitations by studying a hierarchical Gaussian