By employing Monte Carlo random sampling, traditional binary population synthesis (BPS) offers a substantial improvement in efficiency over brute force, grid-based studies. Even so, BPS models typically require a large number of simulation realizations, a computationally expensive endeavor, to generate statistically robust results. Recent advances in statistical methods have led us to revisit the traditional approach to BPS. In this work we describe our publicly available code dart _ board, which combines rapid binary evolution codes, typically used in traditional BPS, with modern Markov chain Monte Carlo methods. dart _ board takes a novel approach that treats the initial binary parameters and the supernova kick vector as model parameters. This formulation has several advantages, including the ability to model either populations of systems or individual binaries, the natural inclusion of observational uncertainties, and the flexible addition of new constraints, which are problematic to include using traditional BPS. After testing our code with mock systems, we demonstrate the flexibility of dart _ board by applying it to three examples: (i) a generic population of high-mass X-ray binaries (HMXBs); (ii) the population of HMXBs in the Large Magellanic Cloud (LMC), in which the spatially resolved star formation history is used as a prior; and (iii) one particular HMXB in the LMC, Swift J0513.4−6547, in which we include observations of the system’s component masses and orbital period. Although this work focuses on HMXBs, dart _ board can be applied to a variety of stellar binaries, including the recent detections by gravitational wave observatories of merging compact object binaries.