  The equation of state (EOS) in the limit of infinite symmetric nuclear matter
exhibits an equilibrium density, $n_0 \approx 0.16 \, \mathrm{fm}^{-3}$, at
which the pressure vanishes and the energy per particle attains its minimum,
$E_0 \approx -16 \, \mathrm{MeV}$. Although not directly measurable, the
nuclear saturation point $(n_0,E_0)$ can be extrapolated by density functional
theory (DFT), providing tight constraints for microscopic interactions derived
from chiral effective field theory (EFT). However, when considering several DFT
predictions for $(n_0,E_0)$ together, a discrepancy between model classes
emerges at high confidence levels that each model prediction's uncertainty
cannot explain. How can we leverage these DFT constraints to rigorously
benchmark nuclear saturation properties of chiral interactions? To address this
question, we present a Bayesian mixture model that combines multiple DFT
predictions for $(n_0,E_0)$ using an efficient conjugate prior approach. The
inferred posterior distribution for the saturation point's mean and covariance
matrix follows a Normal-inverse-Wishart class, resulting in posterior
predictives in the form of correlated, bivariate $t$-distributions. The DFT
uncertainty reports are then used to mix these posteriors using an ordinary
Monte Carlo approach. At the 95% credibility level, we estimate $n_0 \approx
0.157 \pm 0.010 \, \mathrm{fm}^{-3}$ and $E_0 \approx -15.97 \pm 0.40 \,
\mathrm{MeV}$ for the marginal (univariate) $t$-distributions. Combined with
chiral EFT calculations of the pure neutron matter EOS, we obtain bivariate
normal distributions for the nuclear symmetry energy and its slope parameter
evaluated at $n_0$: $S_v \approx 32.0 \pm 1.1 \, \mathrm{MeV}$ and $L\approx
52.6\pm 8.1 \, \mathrm{MeV}$ (95%), respectively. Our Bayesian framework is
publicly available, so practitioners can readily use and extend our results.
