A data-driven model of Earth’s magnetosheath is developed by training a Bayesian recurrent neural network to reproduce Magnetospheric MultiScale (MMS) measurements of the magnetosheath plasma and magnetic field using measurements from the Wind spacecraft upstream of Earth at the first Earth-Sun Lagrange point (L1). This model, called PRIME-SH in reference to its progenitor algorithm PRIME (Probabilistic Regressor for Input to the Magnetosphere Estimation), is shown to predict spacecraft observations of magnetosheath conditions accurately in a statistical sense with a continuous rank probability score (CRPS) of $0.227\sigma$ and more accurately than current analytical models of the magnetosheath. Furthermore, PRIME-SH is shown to reproduce physics not explicitly enforced during training, such as field line draping, the dayside plasma depletion layer, the magnetosheath flow stagnation point, and the Rankine-Hugoniot MHD shock jump conditions. PRIME-SH has the additional benefits of being computationally inexpensive relative to global MHD simulations, being capable of reproducing difficult-to-model physics such as temperature anisotropy, and being capable of reliably estimating its own uncertainty to within $3.5\%$.