Many problems in contemporary astrophysics—from understanding the formation of black holes to untangling the chemical evolution of galaxies—rely on knowledge about binary stars. This, in turn, depends on the discovery and characterization of binary companions for large numbers of different kinds of stars in different chemical and dynamical environments. Current stellar spectroscopic surveys observe hundreds of thousands to millions of stars with (typically) few observational epochs, which allows for binary discovery but makes orbital characterization challenging. We use a custom Monte Carlo sampler (The Joker) to perform discovery and characterization of binary systems through radial velocities, in the regime of sparse, noisy, and poorly sampled multi-epoch data. We use it to generate posterior samplings in Keplerian parameters for 232,495 sources released in APOGEE Data Release 16. Our final catalog contains 19,635 high-confidence close-binary (P ≲ few years, a ≲ few $\mathrm{au}$) systems that show interesting relationships between binary occurrence rate and location in the color-magnitude diagram. We find notable faint companions at high masses (black hole candidates), at low masses (substellar candidates), and at very close separations (mass-transfer candidates). We also use the posterior samplings in a (toy) hierarchical inference to measure the long-period binary-star eccentricity distribution. We release the full set of posterior samplings for the entire parent sample of 232,495 stars. This set of samplings involves no heuristic "discovery" threshold and therefore can be used for myriad statistical purposes, including hierarchical inferences about binary-star populations and subthreshold searches.