Geomagnetic polarity time scales (GPTSs) have been constructed by interpolating between dated marine magnetic anomalies assuming uniformly varying spreading rates. A strategy to obtain an optimal GPTS is to minimize spreading rate fluctuations in many ridge systems; however, this has been possible only for a few spreading centers. We describe here a Monte Carlo sampling method that overcomes this limitation and improves GPTS accuracy by incorporating information on polarity chron durations estimated from astrochronology. The sampling generates a large ensemble of GPTSs that simultaneously agree with radiometric age constraints, minimize the global variation in spreading rates, and fit polarity chron durations estimated by astrochronology. A key feature is the inclusion and propagation of data uncertainties, which weigh how each piece of information affects the resulting time scale. The average of the sampled ensemble gives a reference GPTS, and the variance of the ensemble measures the time scale uncertainty. We apply the method to construct MHTC12, an improved version of the M-sequence GPTS (Late Jurassic-Early Cretaceous, ~160–120 Ma). This GPTS minimizes the variation in spreading rates in a global data set of magnetic lineations from the Western Pacific, North Atlantic, and Indian Ocean NW of Australia, and it also accounts for the duration of five polarity chrons established from astrochronology (CM0r through CM3r). This GPTS can be updated by repeating the Monte Carlo sampling with additional data that may become available in the future.