Determining the timing of paleoearthquakes is a central goal of paleoseismology and serves as an important input to seismic hazard evaluations. Herein, we develop a Bayesian statistical method for refining the ages of strata and earthquakes and calculating the recurrence of earthquakes based on data from paleoseismic excavations. Our work extends previous paleoseismic Bayesian modeling by simultaneously calculating the joint posterior probability density of recurrence intervals, earthquake ages, and layer ages. To estimate this joint posterior density, we employed the Metropolis–Hastings Markov-Chain Monte Carlo (MCMC) sampling method. We used the Wrightwood paleoseismic excavation site as a test bed for our method by comparing the results derived from our simulation to those of previous studies. We found that for stratigraphic order constraints, our results agreed well with previously used methods; however, when peat growth was used to constrain the minimum time that separates layers from one another, a systematic aging of strata in the lower portion of the stratigraphic column resulted. Inspection of our results indicated that, rather than representing a minimum time that separates layers, peat growth instead determined the strata’s ages exactly when a constant peat-growth rate was prescribed. Using this information, we developed a new method of improving layer age, earthquake age, and recurrence interval estimates at sites where peat-bearing layers are present. This method generally produces layer-age estimates with less variance than previous methods and requires that both proximal stratigraphic constraints and trench-wide constraints be obeyed. The benchmarking of the MCMC simulation results shows that simulation provides an effective and efficient way of improving estimates of layer ages, earthquake ages, and recurrence intervals.