The most promising probe into compact binary formation and evolution will be the populations resolved by both gravitational wave (GW) and electromagnetic (EM) observatories. As with any observation, theoretical models must be characterized to gain a fuller understanding of the underlying physics governing the observed population. Investigation of these models in population synthesis studies requires a study of the variance in the observable population from a statistical sample of population realizations. Traditional population synthesis studies have failed to investigate this variance since sizable computational resources are needed to generate a single population. My thesis research has been focused on the development of the Compact Object Synthesis and Monte-Carlo Investigation Code (COSMIC). COSMIC 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 observable compact binary populations using the same computational power as previous population synthesis methods. COSMIC has been used in four studies to date to simulate populations of binary black holes, black holes orbiting luminous stellar companions, and accreting double white dwarfs.