Gravitational-wave (GW) observations of binary neutron star (BNS) mergers can be used to measure luminosity distances and hence, when coupled with estimates for the mergers' host redshifts, infer the Hubble constant H0. These observations are, however, affected by GW measurement noise, uncertainties in host redshifts and peculiar velocities, and are potentially biased by selection effects and the misspecification of the cosmological model or the BNS population. The estimation of H0 from samples of BNS mergers with optical counterparts is tested here by using a phenomenological model for the GW strains that captures both the data-driven event selection and the distance-inclination degeneracy, while being simple enough to facilitate large numbers of simulations. A rigorous Bayesian approach to analyzing the data from such simulated BNS merger samples is shown to yield results that are unbiased, have the appropriate uncertainties, and are robust to model misspecification. Applying such methods to a sample of N≃50 BNS merger events, as LIGO+Virgo could produce in the next ∼5  years, should yield robust and accurate Hubble constant estimates that are precise to a level of ≲2  km s−1 Mpc−1, sufficient to reliably resolve the current tension between local and cosmological measurements of H0.