With the advent of gravitational wave (GW) astronomy, following the first GW detection in September 2015 by the LIGO Scientific Collaboration and the Virgo Collaboration, we now have a means to measure the properties of GW emitters, including merging compact objects such as black holes and neutron stars. The expected number of merger events in a given volume of space and observation time can be modeled as an inhomogeneous Poisson process, with an overall rate and a probability distribution over the systems’ various physical properties. This rate and distribution will provide information on the population of binary star systems, and can be compared with theoretical population models. In this talk, we describe a general Bayesian framework for inferring this rate and distribution, using Markov chain Monte Carlo (MCMC). This method works on noisy data, and relies on a measure of the volume in space to which our detectors are sensitive (which is itself a function of the physical properties). This method relies on a parameterization of the probability distribution of system parameters, which can be restrictive (e.g., a power law) or flexible (e.g., a Gaussian mixture model). With some slight modifications, this method may be used in a variety of terrestrial applications.