MatCalib: A Matlab software package for Bayesian calibration of radiocarbon ages subject to temporal order constraints

Published: 10 September 2021| Version 1 | DOI: 10.17632/rx478cbpm5.1
Contributor:
Shiyong Yu

Description

Radiocarbon ages must be calibrated due to the remarkable fluctuations of the atmospheric radiocarbon level. However, the map from the radiocarbon age domain to the calendar age domain is not one-on-one, providing that the calibration curve is not an injective function. The traditional method only calibrates the radiocarbon age individually, without considering their temporal/stratigraphical ordering. Bayesian radiocarbon age modeling is advantageous over the traditional method in several aspects. First, it can provide a more precise age estimate than the individual calibration by applying some constrains known a priori. Second, it may provide age estimates for an archaeological feature or a geological event that is unable to be dated directly. Third, it represents an adaptable method of statistical inference. According to Bayes’ theorem, the prior information can be formulated mathematically and integrated into the process of inference, which can be easily implemented using the Markov chain Monte Carlo method. Here, a hierarchical Bayesian model with a minimum level of structural complexity is presented. It provides users with a flexible and powerful framework to assemble radiocarbon ages into a sequence along a one-dimensional continuum so that it best reveals their temporal ordering, thereby yielding a more precise timing. The accompanying Matlab software package not only complements the existing MatCal package designed to calibrate radiocarbon ages individually, but also serves as an alternative to the online tools of radiocarbon age calibration such as OxCal and BCal.

Files

Steps to reproduce

Running the software package is quite straightforward. The only file users need to modify is MatCalib_main, which provides a textual interface enabling users to input their data, define their project-specific problem, and run the model to obtain and visualize the results. In the Matlab environment, users need to open the main program MatCalib_main first, supply information about their data in each phase of the sequence including laboratory codes, radiocarbon ages and errors, reservoir ages and errors, depths, and the temporal ordering of the ages (e.g., “ordered”, “unordered”, or “coeval”), and then define the relationship of the phases in terms of “overlapping”, “contiguous”, or “disjoint”, specify the calibration curve (e.g., IntCal13 or IntCal20), give the nearest years, and choose the time scale that the modeled calendar ages are to be rounded and reported. Upon completion of the model run, the empirical probability density function as well as the 68.2 % and 95.4 % HPD regions of the modeled calendar ages of the radiocarbon ages and the phase boundaries are estimated, and the results are saved in files and also displayed in the workspace. Users also can compare any two ages by calculating their posterior difference, and the results are also given in terms of the empirical probability density function as well as the 68.2 % and 95.4 % HPD regions. It is noteworthy that the MCMC method is simulation-based, which requires a large number of iterations to produce reliable (i.e., converged and well mixed) results. Therefore, the model may yield slightly different results at each time of run. Another caveat is that models may be so inconsistent with the data that no possible solutions will be yielded, particularly in case of the presence of extreme outliers. Therefore, users should choose and remove the outliers before running the model to produce valid results. In other cases, the solution space may be very fragmentary, providing that the MCMC solutions will be different with each different starting position. Keeping these properties of the MCMC method in mind, users must check the reproducibility by making multiple runs of the model with different number of iterations, burn-in periods, thinning sizes, and initial values.