MatCalib: a Matlab software package for Bayesian modeling of radiocarbon ages subject to temporal order constraints

1.   Introduction

Radiocarbon (¹⁴C) dating represents a key chronological method that has been widely used for determining the age of geological and archaeological samples of the last ca. 55,000 years ^[1]. The physical basis of this method is that organisms assimilate ¹⁴C, directly or indirectly, from the atmosphere during their life span. However, their death cuts off the linkage to the contemporaneous global carbon cycle, which leaves the remainder ¹⁴C atoms inside their bodies to decay without replenishment with a half-life of ca. 5730 years. Therefore, by the law of radioactive decay, the radiocarbon age of a sample can be determined uniquely by measuring the residual ¹⁴C radioactivity ^[2], providing that the atmospheric ¹⁴C level is assumed to remain constant over time. However, owing to the complex global carbon cycle and the changes in the intensity of the solar and earth's magnetic field ^[3], the atmospheric ¹⁴C level fluctuated remarkably ^[4], exerting a profound impact on the radiocarbon dating method ^[5]. Therefore, the laboratory radiocarbon ages must be calibrated using a calibration curve ^[6], which is a key step for converting radiocarbon ages to calendar ages'.

Due to the non-monotonic nature of the calibration curves, the map from the radiocarbon year to the calendar year is not one to one, thereby leading to non-unique solutions usually expressed in terms of a probability density function with one or several highest posterior density (HPD) regions ^[7,8]. It appears that uncertainty is a nature of the calibrated age, which cannot be eliminated regardless of how small the laboratory error is ^[5]. Sometimes the uncertainty of the calibrated age could be extremely large if the radiocarbon age hits a plateau in the calibration curve ^[9,10]. Many attempts have been made to reduce the uncertainty of radiocarbon age calibration. Bayesian methods represent a promising approach to the solution of this problem ^[11,12,13], which takes advantage of the stratigraphical context or temporal order of the radiocarbon ages as a constraint ^[14]. Specifically, if the radiocarbon ages are assumed to come from a depositional sequence, they must follow a chronological order ^[15,16]. By applying the constraint of temporal order on the calibration, the uncertainty of the calibrated age can be greatly narrowed down ^[8].

Bayesian methods have been increasingly used for modeling a suite of ¹⁴C ages subject to temporal order constraints since 1990s. Several Bayesian frameworks have been developed and in use. Of particular note are the standalone package DateLab ^[17] and the online OxCal ^[11] and BCal ^[18] tools. DateLab provides a simple graphical user interface, while OxCal and BCal allow registered users to access through a worldwide web (www) browser for data entry and workflow setup. However, further analyses of the modeled ages require the random samples from the original Markov Chain Monte Carlo (MCMC) output. For example, if two ages need to be compared, the best way is to calculate the difference between their posterior ages and expresses it in terms of a probability density function with the HPD regions. Also, inexperienced users may need not only a more flexible and user-friendly environment for data entry and definition of their project-specific problems, but also a powerful post-processing tool for analyzing and visualizing the results. Towards this end, a hierarchical Bayesian framework for ¹⁴C age modeling subject to temporal order constraints is described and a Matlab package is presented here. The package is an extension to the previous version ^[19] by adding several new features such as the overlapping boundary relationship that deal with the coexistence and subsequent replacement of one period by another.

2.   Model description

Considering a chronological sequence consisting of M periods, where "period" is defined as a stage of deposition or human occupation characterized by a distinct geological or archaeological feature. Similar to the "phase" in OxCal, each period contains a group of at least one radiocarbon age. For the ith period of the chronological sequence, where $i = 1, 2, \cdots , M$ denoting the number of the lowest (oldest) to the uppermost (most recent) period, Let ${\bf{r}}_{i} = \left[{r}_{i, 1}, {r}_{i, 2}, \cdots , {r}_{i, {N}_{i}}\right]$, $\pmb{\mathsf{σ}}_{i} = \left[{\sigma }_{i, 1}, {\sigma }_{i, 2}, \cdots , {\sigma }_{i, {N}_{i}}\right]$, $\pmb{\mathsf{θ}}_{i} = \left[{\theta }_{i, 1}, {\theta }_{i, 2}, \cdots , {\theta }_{i, {N}_{i}}\right]$, α_i, and β_i denote the N_i ¹⁴C ages, the associated one standard deviation of the N_i ¹⁴C ages, the corresponding calendar ages of the N_i ¹⁴C ages, and the calendar age of the early and late boundary of the ith period, respectively (Figure 1). Note that the division of periods in the sequence is subjective, but mainly depends on the sedimentary facies in geology or tool assemblage in archaeology or the scientific questions to be answered. Here a hierarchical Bayesian framework of ¹⁴C age modeling with a minimum level of structural complexity is described (Figure 2), within which θ_i, α_i, and β_i, where $i = 1, 2, \cdots , M$, are treated as model parameters, which will be inferred using the laboratory radiocarbon ages. Due to lack of knowledge about these parameters, vague (noninformative) priors are preferred to use. For the sake of mathematical expression, the "BP" scale is used for model description, while the algorithms were implemented in both "BP" and "BC/AD" scales. The conversion from the "BP" to the "BCE|CE" scale is BCE|CE = 1950 – BP.

2.1.   Specification of priors

2.1.1.   Calendar age of ¹⁴C ages

For the N_i radiocarbon ages in the ith period, where $i = 1, 2, \cdots , M$, the corresponding calendar ages, say θ_i, are assumed to follow the uniform distribution supported on the interval of $\left[{\alpha }_{i}, {\beta }_{i}\right]$.

where $I\left(\pmb{\mathsf{θ}}_{i}|{\alpha }_{i}, {\beta }_{i}\right)$ is an indicator function describing the temporal order of the calendar ages in the ith period. Here, three types of temporal order are considered: (1) "unordered", which means that the ages are not necessarily to be chronologically ordered in a period. An example is that ages were obtained from different archaeological sites of the same cultural period. As such, the indicator function can be defined as:

(2) "coeval", which means that the ages are assumed to be the same in a period. An example is that ages were obtained from the same depth of a stratigraphical unit or the same archaeological feature of a cultural period and thus their calendar ages must be the same. In this case, the indicator function can be defined as:

and (3) "ordered", which means that the ages follow a chronological order. An example is that ages were obtained from different depths of a stratigraphical unit. Therefore, the indicator function can be expressed as:

where " > " denotes "older than" and "=" denotes "same as".

2.1.2.   Calendar ages of the boundaries

The calendar age of the early and late boundaries of the ith period, say α_i and β_i, where $i = \mathrm{1, 2}, \cdots , M$, is assumed to follow the uniform distribution supported on the interval [A, B], where A and B are two floating parameters defining the beginning and the end of the chronological sequence in calendar years BP, respectively. Therefore, for the first and last periods (i.e., i = 1 or M), the prior probability density function of the calendar ages of the early and late boundaries, say α_i and β_i, respectively, can be written as:

where $I\left({\alpha }_{i}, {\beta }_{i}|A, B\right)$ is an indicator function describing the temporal order of the calendar ages of the early and late boundaries. It can be simply defined as:

However, for the internal periods, the prior probability density function of the calendar ages of the early and late boundaries should be treated differently according to the relationship with the neighboring periods. Without loss of generality, for the ith internal period, where $i = 2, \cdots , M-1$, the prior probability density function of the early and late boundaries, say α_i and β_i, respectively, can be written as:

where $I\left({\alpha }_{i}, {\beta }_{i}|{\alpha }_{i-1}, {\beta }_{i-1}, {\alpha }_{1}, {\beta }_{M}\right)$ is an indicator function describing the relationship of two neighboring periods. This mathematical framework is essentially the same as what Nicholls and Jones ^[14] proposed. Here three cases are considered: (1) "overlapping", which means the imbrication of two neighboring periods. An example is the initial coexistence of two archaeological cultures and subsequent replacement of one culture by another. Therefore, the indicator function can be expressed as:

(2) "contiguous", which means a sudden transition from one period to another. An example is the rapid change of the depositional environment or the archaeological culture. As such, the indicator function can be defined as:

and (3) "disjoint", which means that there is an interruption between two periods. An example is the existence of a depositional hiatus between two stratigraphical units or a gap of human occupation between two archaeological periods. As such, the indicator function can be written as:

where " > " denotes "older than" and "=" denotes "same as".

2.2.   Specification of the likelihood

The normal distribution is usually used for modeling ¹⁴C ages from the calibration curve ^[20,21]. However, the traditional normal distribution is less robust to outlying data. To better accommodate outliers in the ¹⁴C ages, we adopt the variance-multiplier normal distribution ^[22]. Let µ and ϵ be the mean and one standard deviation of a radiocarbon calibration curve, respectively. For the jth calendar age in the ith period, say ${\theta }_{i, j}$, where $i = 1, 2, \cdots , M$ and $j = 1, 2, \cdots , {N}_{i}$, the corresponding true ¹⁴C age is assumed to be normally distributed around the calibration curve mean, $\pmb{\mathsf{μ}}\left({\theta }_{i, j}\right)$, confined by an additive variance of the laboratory ¹⁴C age and the calibration curve multiplied by a positive number. Hence, the true ¹⁴C age can be modeled as:

where ${{\omega }_{i, j}^{2} = \sigma }_{i, j}^{2}+\pmb{\mathsf{ϵ}}^{2}\left({\theta }_{i, j}\right)$ is the sum of the laboratory ¹⁴C age variance (i.e., ${\sigma }_{i, j}^{2}$) and the calibration curve variance (i.e., $\pmb{\mathsf{ϵ}}^{2}\left({\theta }_{i, j}\right)$) and δ is a positive number. Let δ be a random variable with an improper prior probability density function defined as:

where $\delta \in \left(0, +\infty \right)$. As such, for the ith period, where $i = 1, 2, \cdots , M$, the likelihood function can be expressed as:

An integrated likelihood can be obtained by integrating out hyperparameter δ (Supplementary), yielding:

where Γ (·) denotes the gamma function.

2.3.   Combining the priors and likelihood to estimate the posterior

Making use of Bayes' theorem and the integrated likelihood function (Eq. 14), for the ith period, where $i = 1, 2, \cdots , M$, the joint posterior probability density function of parameters θ_i, α_i, and β_i, can be expressed as:

Note that the prior indicator functions of the calendar ages and the boundaries have different expressions according to the ordering of the ages and the relationship of the boundaries as given above.

3.   Algorithm

The structure of the above equation (Eq. 15) is so complex that we are unable to derive an analytical expression of the marginal posterior probability distribution function for each parameter. Therefore, this model framework can be implemented using the MCMC method. The Metropolis-approximated Gibbs sampler is used to generate a Markov Chain for the parameters one at a time, and a Metropolis-Hastings step is used to update each parameter ^[23]. Specifically, for the ith period of the sequence, where $i = 1, 2, \cdots , M$, we choose to update the parameters in the order θ_i → α_i → β_i. Given the current state of the calendar ages of the ¹⁴C ages, say θ_i, and the early and late boundaries, say α_i and β_i, respectively, we first update θ_i using a uniform random walk process on the interval [α_i, β_i] such that:

where u denotes a vector of N_i independent random numbers drawn from a uniform distribution on [–1, 1], and ρ is the step size of the random walk. Then, we update the calendar age of the early boundary of the ith period of the sequence, say α_i, where $i = 1, 2, \cdots , M$, using the triangular distribution. For the first period, the calendar age of the early boundary, say α₁, can be updated according to:

where T (·) denotes the triangular distribution with three parameters describing the lower bound, the peak, and the upper bound, respectively. While for other periods (i.e., i > 1), special treatments are required, depending on the relationship of two neighboring periods. Three cases are considered accordingly. For the "overlapping" boundaries, the calendar age of the early boundary can be updated according to:

while for the "contiguous" boundaries, the calendar age of the early boundary can be updated according to:

and for the "disjoint" boundary, the calendar age of the early boundary can be updated according to:

Finally, we update the calendar age of the late boundary of the ith period, say β_i, where $i = 1, 2, \cdots , M$, using the triangular distribution too. Similarly, for the first M − 1 periods, three cases are considered. For the "overlapping" boundaries, the calendar age of the late boundary can be updated according to:

while for the "contiguous" boundaries, the calendar age of the late boundary, say β_i, can be updated by simply setting:

and for the "disjoint" boundaries, the calendar age of the late boundary can be updated according to:

For the last period (i.e., i = M), the calendar of the late boundary, say β_M, can be updated according to:

The acceptance or rejection of the proposed move of the chains is determined using the Metropolis-Hastings algorithm ^[24]. The state-of-the-art IntCal20 ^[25], SHCal20 ^[26], and Marine20 ^[27] calibration datasets are also provided for modeling of the ¹⁴C ages. For the sake of computing, all of the laboratory and the calibration curve radiocarbon ages are converted to the F¹⁴C space during the modeling process ^[20]. The convergence of the chains is monitored using the method proposed by Gelman ^[28]. Several chains are run parallelly with different starting points for each set of simulations, and the ratio of the between-sequence to with-sequence variance in terms of the potential scale reduction factor, $\widehat{R}$ is calculated. If $\widehat{R}$ is below the critical value of 1.1, the chains can be regarded as convergent. Once the chains for a model parameter converge, they are mixed to allow the estimate of the empirical probability density function that mimics the posterior probability distribution of this parameter. We implemented the above procedure in the Matlab^® environment, yielding a software package named as MatCalib v2.0. The code is available on Mendeley Data ^[29] and also accessible in the Supplementary.

4.   User guide

Running the software package is quite straightforward. A tab delimited spread sheet template is provided, which allow users to enter their data, including laboratory codes, depths, ¹⁴C ages and errors, reservoir ages and errors if any, specify the calibration curve (e.g., IntCal13 or IntCal20) for each age, and define the temporal ordering of the ages (e.g., "ordered", "unordered", or "coeval") in each period. The only file users need to modify is MatCalib_main.m, which provides a textual interface enabling users to read their data, define their project-specific problem (e.g., boundary relationships), specify parameters for MCMC simulations, and run the model to obtain, analyze, and visualize the results. In the Matlab environment, users need to open the main program MatCalib_main, give the nearest years that the modeled calendar ages are to be rounded to, choose the time scale in which the modeled calendar ages are to be reported, and specify the parameters for Monte Carlo simulations such as the number of chains to be run for convergence diagnoses as well as the length, burn-in period, and thinning interval of the chains. Then, MCMC simulations are conducted using the function age_modeling and the convergence of the chains is tested using the function convergence.

Upon completion of the MCMC simulations, the empirical probability density function, as well as the 68.2% and 95.4% HPD regions of the modeled calendar ages of the ¹⁴C ages and the early and late boundaries of each period, are estimated using the method of Lougheed and Obrochta ^[30], and the results are saved automatically to files for subsequent analyses. Users can plot the results either against sample numbers using the function plot_ages, or against depths using the function plot_age_depth if depth information is provided. Users also can compare any two ages by calculating their posterior difference and calculate the pooled mean of several ages using the functions age_difference and pooled_mean, respectively. The results are also given in terms of the empirical probability density function as well as the 68.2% and 95.4% HPD regions, and they can be plotted using the functions plot_difference and plot_pooled_mean, respectively.

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 a run. Another caveat is that models may be so inconsistent with the data that no possible solutions will be yielded, particularly in the 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 a different number of iterations, burn-in period, thinning size, and initial value (automatically reset for each run using the function initialization).

5.   Illustrative examples

Some of the functionalities of this software package are demonstrated using two datasets of ¹⁴C ages. One is obtained from the Neolithic archaeological sites on the Shandong Peninsula in North China, which is used to determine the timing of the Longshan and Yueshi cultures in this area. The other is retrieved from a sedimentary core on the southeastern Swedish Baltic coast, which is used to establish an age-depth relationship and to quantify the sedimentary hiatus induced by an erosional event during the transition from a lacustrine to a lagoonal environment. In both cases, a set of simulations with three chains are run for 20,000 iterations for each parameter. Every 20^th iteration is kept after a burn-in period of 1000 iterations to reduce autocorrelation.

5.1.   Timing of the Longshan and Yueshi cultures on the Shandong Peninsula

The dataset used in this study consists of nine ¹⁴C ages from sites belonging to the Longshan cultural period and 10 ages belonging to the Yueshi cultural period ^[31]. These ages are put together into a chronological sequence with two periods corresponding to the Longshan (Period 1) and Yueshi (Period 2), respectively. The latest Longshan site was dated to 3570 ± 80 ¹⁴C yr BP and the earliest Yueshi site was dated to 3760 ± 145 ¹⁴C yr BP, implying that these two cultures may have coexisted for a while ^[32]. Therefore, overlapping boundaries are used in the age modeling. As the ¹⁴C ages in each cultural period come from individual sites and their stratigraphical relationship is unknown, the temporal order of the ages is set to "unordered". The probability distributions of the modeled calendar ages are presented in Figure 3. The 95.4% HPD regions of calendar age of the boundaries are listed in Table 1 for the sake of comparison with those obtained from OxCal.

The Longshan culture, colloquially referred to as the black pottery culture, is a late-Neolithic culture in the middle and lower Yellow River areas of northern China from about 3000 to 1900 BCE. The culture flourished dramatically and expanded to the Shandong Peninsula at about 2885–2750 BCE as defined by the 95.4% HPD region (Table 1). It then decreased and vanished around 2005–1870 BC (95.4% HPD region). The Yueshi culture rose locally at about 2490–2335 BC (95.4% HPD region). It then replaced the Longshan culture and gradually evolved into the Bronze Age at about 1480–1300 BC (95.4% HPD region). These timings lie in the 95.4% age ranges obtained from OxCal (Table 1), suggesting that MatCalib can yield age estimates much firmer than OxCal. As defined by the age difference between the early boundary of Period 2 (i.e., Yueshi period) and the late boundary of Period 1 (i.e., Longshan period) (Figure 3), these two cultures may have coexisted for about 350–590 years (95.4% HPD region) on the Shandong Peninsula (Figure 4). Again, OxCal yields a relatively large overlap (280–740 years, 95% HPD region) of these two cultures.

5.2.   Quantifying the sedimentary hiatus on the southeastern Swedish Baltic coast

Previous work about the ¹⁴C dating of a sediment core from the Hunnemara ancient lake (14°53' E, 56°10' N) on the southeastern Swedish Baltic coast revealed a depositional hiatus across the transition from the lacustrine to the lagoonal facies ^[33], which may have resulted from erosion during the rapid sea-level rise about 7600 years ago ^[34]. The lake is developed within a closed basin with an outflow threshold around 3.0 m as1. It has been drained and converted to agricultural land and a meadow by local farmers in the 19th century and is now partly used for a garbage dump. An 8.5-m-long core was taken near the center of the lake basin. A total of 10 ¹⁴C ages on bulk sediments were obtained, which are incorporated into a chronological sequence consisting of two periods representing the lacustrine (Period 1) and lagoonal (Period 2) environments, respectively. Period 1 comprises three ¹⁴C ages, whereas Period 2 contains seven ¹⁴C ages subject to a local reservoir age offset of −108 ± 24 years. Therefore, the Marine20 calibration curve is applied on these ages ^[27]. As the ages are obtained from different depths along the core, the law of superposition was strictly enforced by using a temporal order constraint (e.g., "ordered") on the ages in both periods. The hiatus lies at 575 cm, which is marked by a non-depositional surface separating the two periods. Therefore, a "disjoint" relationship of the internal boundaries is applied.

The modeled calendar ages of the radiocarbon ages are plotted against their depths (Figure 5). The age-depth relationship reveals variable sedimentation rates within the lake basin. The results suggest that the lacustrine environment (Period 1) prevailed from 11,745–11,545 BP (95.4% HPD region) to 8530–8145 BP (95.4% HPD region), while the lagoonal environment (Period 2) occurred in the basin at about 7410–7040 BP (95.4% HPD region) as a result of rapid sea-level rise ^[34] and terminated at about 3090–2925 BP (95.4% HPD region). As defined by the posterior age difference of the late boundary of Period 1 and the early boundary of Period 2, the sedimentary hiatus is estimated to be 870–1355 years (95.4% HPD region) (Figure 6).

6.   Discussion

Bayesian methods are widely used for modeling ¹⁴C ages. It represents a flexible yet robust approach to integrating prior knowledge about the ages and their temporal or stratigraphical context, thereby yielding a more precise timing of the geological or archaeological events than the traditional methods ^[8]. Although some key components of the Bayesian methods for ¹⁴C age modeling have already been well explained ^[8,14,21], a formal treatment of several other outstanding issues (e.g., the relationship of neighboring boundaries) is lacking. Here a complete and detailed description of the analytical framework is presented along with a Matlab package for software implementation. Therefore, it can be used not only as an introductory guild to, but also a tutorial on Bayesian ¹⁴C age modeling.

There exist several software packages designed to perform Bayesian ¹⁴C age modeling ^[11,17,18]. Most of which shares some common features. For example, the law of superposition is the backbone of all Bayesian methods, which has been strictly enforced to observe the constraint of temporal order. Yet the model presented here differs from these tools in several aspects. First is the treatment of the relationship between two neighboring periods. As illustrated by the first example, coexistence of two cultures appears to be a common phenomenon in archaeology ^[32]. Similar to OxCal, an overlapping relationship is considered in this model, while this functionality was not implemented in other tools. Second is the application of different calibration curves on ¹⁴C ages of different types. This is particularly necessary if ages were obtained from different depositional environments (e.g., marine versus terrestrial). As illustrated by the second example, ages in the first period are obtained from a terrestrial setting, whereas those in the second period are retrieved from a marine environment subject to the reservoir effect ^[35]. As such, a local reservoir age and the marine calibration curve should be applied on these ages. Other minor differences of this model from the existing tools include the use of the non-Gaussian likelihood function to accommodate outlying data, the implementation of convergence diagnosis of the Markov Chains, and the employment of a triangular distribution to update the boundaries during the simulations.

The accompanying Matlab package (i.e., MatCalib 2.0) contains a wide range of functionality for Bayesian statistical analyses of ¹⁴C ages, including the inference of the calendar age of the individual ¹⁴C ages and the boundary of the periods. Random samples of these parameters generated from the MCMC simulations are saved automatically to a data file for future use without running the model again. Based on these data, the empirical probability density functions of these parameters are built up first and then the HPD regions and the median-probability age are calculated. Further analyses such as pooled mean of ages, timing, and duration of an event can be conducted and visualized using the high-quality graphical output (e.g., Figures 3–6). In future, the software will be improved to enable the modeling of multi-source ages, for example, calendar ages with errors and/or tie points for events without errors. For these age types, a calibration curve is no longer needed. If depth information is available, an age-depth relationship may be established through polynomial regression (e.g., CLAM ^[36]), Brownian bridge (e.g., OxCal ^[37]), or interpolation (e.g., Undatable ^[38]) on the stored random samples at known depths, instead of employing a stochastic depositional model such as Bchron ^[39] and Bacon ^[40].

7.   Conclusions

The traditional method of ¹⁴C age calibration maps individual ¹⁴C ages to the calibration curve to obtain the calendar age. However, owing to the large fluctuations of the atmospheric ¹⁴C level, the resulting calibrated age is not a point estimate; rather, it is usually treated as a random variable expressed in terms of a probability distribution with one or several HPD regions. Most importantly, the probability distribution of the calibrated age could be extremely broad if the radiocarbon age hits a radiocarbon plateau in the calibration curve, thereby yielding several large HPD regions. The large uncertainty in the calibrated age hinders the precise timing of geological and/or archaeological events. Conceptually different from the traditional method, Bayesian model provides such an analytical framework that integrate prior knowledge about the ages and their temporal order as a constraint, which in turn can significantly narrow down the uncertainty of the calibrated ages.

The hierarchical Bayesian model presented in this paper provides a flexible approach to the precise timing of archaeological and/or geological events by making use of a suite of ¹⁴C ages and their stratigraphical context (temporal order). Implemented in the Matlab environment, the accompanying software package is not only an extension to the existing MatCal package that is only able to calibrate individual ¹⁴C ages ^[30], but also an alternative or complement to the online calibration tools such as BCal ^[18] and OxCal ^[11]. Also, this model can be used as a building block for age-depth modeling. If the depth information of the ¹⁴C ages is available, a probabilistic expression of the age-depth relationship can be readily obtained through either bootstrap regression or Brownian bridge transformation.

Acknowledgments

This work was supported by the National Science Foundation of China (grant number 41971102). My gratitude is due to the three anonymous reviewers for their insightful comments and constructive suggestions that greatly improved the manuscript.

Conflict of interest

The author declares no conflict of interest.

Supplementary

1. Derivation of the integrated likelihood function.

2. Matlab code for implementing Bayesian radiocarbon age modeling.