The ubiquitous multiple velocity components seen in turbulent molecular clouds pose a challenge to conventional spectral line analysis techniques. To this end we have developed a novel Bayesian methodology to determine the number of velocity components and estimate their model parameters in a statistically rigorous and automated way. We have written an open-source software implementation for community use, NestFit. The method selects models using Bayes factors computed with Nested Sampling Monte Carlo. Application to the GBT GAS synthetic test suite shows superior retrieval accuracy compared to leading automated methods. We present results from an analysis of the GBT KEYSTONE M17 and Mon R2 fields. The results agree well with previous analyses for spectra at high-SNR, but the new techniques produce substantial improvement at low-SNR due to (i) the numerical robustness of Nested Sampling and (ii) the conditioning of the data on empirically motivated priors. The Bayesian inference methods presented in this work enable reproducible analyses of large observational and simulated datasets in a flexible and automated fashion without being restricted to the specific model and training dataset used, as is the case for Neural Network based approaches. Further, the Bayesian analysis of spectra in the low-SNR regime enabled by NestFit represents a powerful new tool for studying the molecular ISM at low densities, such as: the environments of dense cores, large-scale converging flows, and magnetically sub-critical gas.