In this paper, we present a variational Bayesian (VB) approach to computing the interval estimates for nonhomogeneous Poisson process (NHPP) software reliability models. This approach is an approximate method that can produce analytically tractable posterior distributions. We present simple iterative algorithms to compute the approximate posterior distributions for the parameters of the gamma-type NHPP-based software reliability model using either individual failure time data or grouped data. In numerical examples, the accuracy of this VB approach is compared with the interval estimates based on conventional Bayesian approaches, i.e., Laplace approximation, Markov chain Monte Carlo (MCMC) method, and numerical integration. The proposed VB approach provides almost the same accuracy as MCMC, while its computational burden is much lower.