Over the past twenty years, analyzing the abundance of the isotope chlorine-36 ( 36 Cl) has emerged as a popular tool for geologic dating. In particular, it has been observed that 36 Cl measurements along a fault plane can be used to study the timings of past ground displacements during earthquakes, which in turn can be used to improve existing seismic hazard assessment. This approach requires accurate simulations of 36 Cl accumulation for a set of fault-scarp rock samples, which are progressively exhumed during earthquakes, in order to infer displacement histories from 36 Cl measurements. While the physical models underlying such simulations have continuously been improved, the inverse problem of recovering displacement histories from 36 Cl measurements is still mostly solved on an ad-hoc basis. The current work resolves this situation by providing a MATLAB implementation of a fast, automatic, and flexible Bayesian Markov-chain Monte Carlo algorithm for the inverse problem, and provides a validation of the 36 Cl approach to inference of earthquakes from the demise of the Last Glacial Maximum until present. To demonstrate its performance, we apply our algorithm to a synthetic case to verify identifiability, and to the Fiamignano and Frattura faults in the Italian Apennines in order to infer their earthquake displacement histories and to provide seismic hazard assessments. The results suggest high variability in slip rates for both faults, and large displacements on the Fiamignano fault at times when the Colosseum and other ancient buildings in Rome were damaged.