Making a prediction typically involves dealing with uncertainties. The application of uncertainty analysis to buildings and HVAC (heating, ventilation and air conditioning) systems, however, remains limited. Most existing studies concentrate on the parameter uncertainty and parametric variability in building simulations for the design stage, and rely on Monte Carlo experiments to quantify this uncertainty. This dissertation aims to develop a rapid and direct method that is capable of quantifying uncertainty when predicting building cooling and heating consumption in the operation stage, while simultaneously capturing all sources of uncertainty and applying these to actual system operations. Gaussian Process regression, a Bayesian modeling method, is proposed for this purpose. The primary advantage of Gaussian Process regression is that it directly outputs a probability distribution that explicitly expresses prediction uncertainty. The predictive distribution covers uncertainty sources arising not only from parameter uncertainty and parametric variability, but also from modeling inadequacy and residual variability. By assuming a Gaussian input distribution and using Gaussian kernels, Gaussian Process regression takes parameter uncertainty and parametric variability into consideration without using the Monte Carlo method. This dissertation makes three main contributions. First, based on the observations from commissioning projects for approximately twenty campus buildings, some of the important uncertainties and typical problems in variable air volume system (VAV) operations are identified. Second, Gaussian Process regression is used to predict building cooling and heating consumption and to evaluate the impact of parametric variability of system control related variables. Third, a method for automated fault detection that uses Gaussian Process regression to model baselines is developed. By using the uncertainty outputs from the Gaussian Process regression together with Bayes classifiers and probabilistic graphical models, the proposed method can detect whether system performance is normal or faulty at the system component level or the whole building level with a high degree of accuracy.