High-fidelity radio interferometric data calibration that minimizes spurious spectral structure in the calibrated data is essential in astrophysical applications, such as 21 cm cosmology, which rely on knowledge of the relative spectral smoothness of distinct astrophysical emission components to extract the signal of interest. Existing approaches to radio interferometric calibration have been shown to impart spurious spectral structure to the calibrated data if the sky model used to calibrate the data is incomplete. In this paper, we introduce BayesCal: a novel solution to the sky-model incompleteness problem in interferometric calibration, designed to enable high-fidelity data calibration. The BayesCal data model supplements the a priori known component of the forward model of the sky with a statistical model for the missing and uncertain flux contribution to the data, constrained by a prior on the power in the model. We demonstrate how the parameters of this model can be marginalized out analytically, reducing the dimensionality of the parameter space to be sampled from and allowing one to sample directly from the posterior probability distribution of the calibration parameters. Additionally, we show how physically motivated priors derived from theoretical and measurement-based constraints on the spectral smoothness of the instrumental gains can be used to constrain the calibration solutions. In a companion paper, we apply this algorithm to simulated observations with a HERA-like array and demonstrate that it enables up to four orders of magnitude suppression of power in spurious spectral fluctuations relative to standard calibration approaches.