Bias and correction due to Monte Carlo integrals in hierarchical inference of gravitational wave populations

With the ever-increasing catalogs of gravitational wave observations, interest is shifting away from individual sources to population level constraints. However, the Bayesian analysis of a population requires an empirical approximation to the hierarchical likelihood, which carries some intrinsic uncertainty. For some populations, this uncertainty is small, whereas for others, the uncertainty can become quite large. Here, we show that such procedures give generically biased constraints on the population, and we discuss the mathematical theory behind this issue. We derive a correction which removes the leading-order biasing term, and provide recommendations to population analysts for handling the likelihood uncertainty in the future.