The most promising probe into compact binary formation and evolution will be the populations observed by both gravitational wave (GW) and electromagnetic (EM) telescopes. As with any observation, biases must be characterized to gain a deeper understanding of the underlying physics governing the observed population. Investigation of these biases in population synthesis studies requires a study of the variance in the observable population from a statistical sample of population realizations. In this talk, I will introduce a population synthesis code, COSMIC, which provides a key necessity in understanding the compact binary populations observable by current and future GW and EM observatories: the ability to generate a statistical sample of synthetic compact binary populations using the same computational power as previous population synthesis methods. I will also show the results from three recent studies which used COSMIC to show (1) LISA will help to constrain binary black hole formation channels, (2) LISA and Gaia together can characterize hundreds of accreting double white dwarfs, and (3) Gaia could discover thousands of black holes by the end of its mission.