  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 with mass ratios $q\leq2$, and spin magnitudes $\chi_1,\chi_2\leq0.8$.
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 (improvement is also reported in the
extrapolated region at high mass ratios and spins). In addition, we also model
the remnant spin and kick directions. Being trained directly on precessing
simulations, our fits are free from ambiguities regarding the initial frequency
at which precessing quantities are defined. We also construct a model for
remnant properties of aligned-spin systems with mass ratios $q\leq8$, and spin
magnitudes $\chi_1,\chi_2\leq0.8$. 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.
