  In the absence of numerous gravitational-wave detections with confirmed
electromagnetic counterparts, the "dark siren" method has emerged as a leading
technique of gravitational-wave cosmology. The method allows redshift
information of such events to be inferred statistically from a catalogue of
potential host galaxies. Due to selection effects, dark siren analyses
necessarily depend on the mass distribution of compact objects and the
evolution of their merger rate with redshift. Informative priors on these
quantities will impact the inferred posterior constraints on the Hubble
constant ($H_0$). It is thus crucial to vary these unknown distributions during
an $H_0$ inference. This was not possible in earlier analyses due to the high
computational cost, restricting them to either excluding galaxy catalogue
information, or fixing the gravitational-wave population mass distribution and
risking introducing bias to the $H_0$ measurement. This paper introduces a
significantly enhanced version of the Python package GWCOSMO, which allows
joint estimation of cosmological and compact binary population parameters. This
thereby ensures the analysis is now robust to a major source of potential bias.
The gravitational-wave events from the Third Gravitational-Wave Transient
Catalogue are reanalysed with the GLADE+ galaxy catalogue, and an updated, more
reliable measurement of $H_0=69^{+12}_{-7}$ km s$^{-1}$ Mpc$^{-1}$ is found
(maximum a posteriori probability and 68% highest density interval). This
improved method will enable cosmological analyses with future
gravitational-wave detections to make full use of the information available
(both from galaxy catalogues and the compact binary population itself), leading
to promising new independent bounds on the Hubble constant.
