A novel method is proposed to infer Bayesian predictions of computationally expensive models. The method is based on the construction of quadrature rules, which are well-suited for approximating the weighted integrals occurring in Bayesian prediction. The novel idea is to construct a sequence of nested quadrature rules with positive weights that converge to a quadrature rule that is weighted with respect to the posterior. The quadrature rules are constructed using a proposal distribution that is determined by means of nearest neighbor interpolation of all available evaluations of the posterior. It is demonstrated both theoretically and numerically that this approach yields accurate estimates of the integrals involved in Bayesian prediction. The applicability of the approach for a fluid dynamics test case is demonstrated by inferring accurate predictions of the transonic flow over the RAE2822 airfoil with a small number of model evaluations. Here, the closure coefficients of the Spalart-Allmaras turbulence model are considered to be uncertain and are calibrated using wind tunnel measurements.