A predictive theory of nuclei and nuclear matter requires full, rigorous uncertainty quantification in nuclear physics calculations. This requires that we quantify not only experimental and numerical errors, but also the errors in our theoretical models. Bayesian Model Mixing (BMM) is emerging as a novel approach to reach this aim. BMM allows us to combine multiple physics models using weights that vary point-to-point across the input space, so that each model is leveraged in its region of applicability. This technique brings much needed uncertainty quantification to any chosen physics application, including the dense matter equation of state (EOS), which itself has seen breakthrough developments in recent years. These include astrophysical observations from neutron star mergers, improvements on microscopic many-body calculations and their uncertainties, heavy-ion collision data, and novel truncation error estimates for chiral effective field theory (EFT). Even with these achievements, there remains no microscopically motivated theory that can predict the dense matter EOS from saturation density up to deconfined quark matter; however, chiral EFT, the functional renormalization group (FRG), and perturbative QCD (pQCD) describe their respective local density regions very well. Using BMM software I developed for the BAND collaboration’s cyberinfrastructure framework, I will introduce our multivariate BMM method and present an application to an EFT-based toy model. I will then discuss multivariate BMM applied to chiral EFT and pQCD, show results for the pressure and speed of sound of symmetric matter, and compare to FRG and heavy ion collision data.