We introduce Bayesian Expansion (BE), an approximate numerical technique for passage time distribution analysis in queueing networks. BE uses a class of Bayesian networks to approximate the exact joint probability density of the model by a product of conditional marginal probabilities that scales efficiently with the model size. We show that this naturally leads to decomposing a queueing network into a set of Markov processes that jointly approximate the dynamics of the model and from which passage times are easily computed. Approximation accuracy of BE depends on the specific Bayesian network used to decompose the joint probability density. Hence, we propose a selection algorithm based on the KullbackLeibler divergence to search for the Bayesian network that provides the most accurate results. Random models and case studies of increasing complexity show the significant accuracy gain of distribution estimates returned by BE compared to Markov and Chebyshev inequalities that are fre...