We introduce a new framework for quantifying correlated uncertainties of the infinite-matter equation of state derived from chiral effective field theory (χEFT). Bayesian machine learning via Gaussian processes with physics-based hyperparameters allows us to efficiently quantify and propagate theoretical uncertainties of the equation of state, such as χEFT truncation errors, to derived quantities. We apply this framework to state-of-the-art many-body perturbation theory calculations with nucleon-nucleon and three-nucleon interactions up to fourth order in the χEFT expansion. This produces the first statistically robust uncertainty estimates for key quantities of neutron stars. We give results up to twice nuclear saturation density for the energy per particle, pressure, and speed of sound of neutron matter, as well as for the nuclear symmetry energy and its derivative. At nuclear saturation density the predicted symmetry energy and its slope are consistent with experimental constraints.