  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, $H_0$. 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 mis-specification of the cosmological model or the BNS
population. The estimation of $H_0$ 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 mis-specification.
Applying such methods to a sample of $N \simeq 50$ BNS merger events, as
LIGO+Virgo could produce in the next $\sim 5$ years, should yield robust and
accurate Hubble constant estimates that are precise to a level of $\sim 2$ km
s$^{-1}$ Mpc$^{-1}$, sufficient to reliably resolve the current tension between
local and cosmological measurements of $H_0$.
