Geochronology is essential for understanding Earth’s history. The availability of precise and accurate isotopic data is increasing; hence it is crucial to develop transparent and accessible data reduction techniques and tools to transform raw mass spectrometry data into robust chronological data. Here we present a Monte Carlo sampling approach to fully propagate uncertainties from linear regressions for isochron dating. Our new approach makes no prior assumption about the causes of variability in the derived chronological results and propagates uncertainties from both experimental measurements (analytical uncertainties) and underlying assumptions (model uncertainties) into the final age determination. Using synthetic examples, we find that although the estimates of the slope and y-intercept (hence age and initial isotopic ratios) are comparable between the Monte Carlo method and the benchmark “Isoplot” algorithm, uncertainties from the later could be underestimated by up to 60%, which are likely due to an incomplete propagation of model uncertainties. An additional advantage of the new method is its ability to integrate with geological information to yield refined chronological constraints. The new method presented here is specifically designed to fully propagate errors in geochronological applications involves linear regressions such as Rb-Sr, Sm-Nd, Re-Os, Pt-Os, Lu-Hf, U-Pb (with discordant points), Pb-Pb and Ar-Ar.