  Multi-messenger observations of binary neutron star mergers offer a promising
path towards resolution of the Hubble constant ($H_0$) tension, provided their
constraints are shown to be free from systematics such as the Malmquist bias.
In the traditional Bayesian framework, accounting for selection effects in the
likelihood requires calculation of the expected number (or fraction) of
detections as a function of the parameters describing the population and
cosmology; a potentially costly and/or inaccurate process. This calculation
can, however, be bypassed completely by performing the inference in a framework
in which the likelihood is never explicitly calculated, but instead fit using
forward simulations of the data, which naturally include the selection. This is
Likelihood-Free Inference (LFI). Here, we use density-estimation LFI, coupled
to neural-network-based data compression, to infer $H_0$ from mock catalogues
of binary neutron star mergers, given noisy redshift, distance and peculiar
velocity estimates for each object. We demonstrate that LFI yields
statistically unbiased estimates of $H_0$ in the presence of selection effects,
with precision matching that of sampling the full Bayesian hierarchical model.
Marginalizing over the bias increases the $H_0$ uncertainty by only $6\%$ for
training sets consisting of $O(10^4)$ populations. The resulting LFI framework
is applicable to population-level inference problems with selection effects
across astrophysics.
