We present accurate fits for the remnant properties of generically precessing binary black holes, trained on large banks of numerical-relativity simulations. We use Gaussian process regression to interpolate the remnant mass, spin, and recoil velocity in the 7-dimensional parameter space of precessing black-hole binaries. For precessing systems, our errors in estimating the remnant mass, spin magnitude, and kick magnitude are lower than those of existing fitting formulae by at least an order of magnitude. In addition, we also model the remnant spin and kick directions. Improvement is also reported for aligned-spin systems. Being trained directly on precessing simulations, our fits are free from ambiguities regarding the initial frequency at which precessing quantities are defined. As a byproduct, we also provide error estimates for all fitted quantities, which can be consistently incorporated into current and future gravitational-wave parameter-estimation analyses. Our model(s) are made publicly available through a fast and easy-to-use Python module called surfinBH.