Analyzing Statistical Age Models to Determine the Equivalent Dose and Burial Age Using a Markov Chain Monte Carlo Method

Galbraith et al. (1999) recommended using the CAM to provide a representative estimate of the characteristic dose for samples that have not been affected by partial bleaching or sediment mixing. The benefit of the CAM is that it can take any over-dispersion (OD) into account when estimating the weighted mean dose and its standard error (Galbraith and Roberts, 2012). OD provides an estimate of the degree of dispersion that represents the relative spread in the D_e distribution that remains after allowing for measurement uncertainties. Routinely, OD is calculated as a descriptive statistic to reveal sources of variability in D_e (e.g., Thomsen et al., 2005; Jacobs et al., 2006). The OD of D_e distributions in well-bleached analogues is commonly used as an input for statistical age models when analyzing measured sedimentary samples (Galbraith and Roberts, 2012). The likelihood function of the CAM can be described as follows: (2.1) p(Dj|θ)=CAM(xj,yj|μ,σ)=12π(xj2+σ2)e−(yj−u)22(xj2+σ2) p\left( {{D}_{j}}|\theta \right)=CAM\left( {{x}_{j}},{{y}_{j}}|\mu ,\sigma \right)=\frac{1}{\sqrt{2\pi \left( x_{j}^{2}+{{\sigma }^{2}} \right)}}{{e}^{-\frac{{{\left( {{y}_{j}}-u \right)}^{2}}}{2\left( x_{j}^{2}+{{\sigma }^{2}} \right)}}} where y_j and x_j represent for the jth logged D_e value and its relative standard error (RSE), respectively. p(D_j|θ) represents the probability of data D_j given parameter θ. μ and σ represent the parameters of interest (i.e., the logged characteristic D_e and the OD, respectively).

Optimal parameters can be obtained using either an iterative numeric procedure or by direct numerical maximization of the joint-likelihood determined from paired observations D={x_j, y_j | j=1,2,3,...,n) (Galbraith et al., 1999). Most free and commercial optimization environments (such as R, Matlab, etc.) provide well-designed optimization routines for optimization usage. The standard errors of parameters are conventionally obtained by inverting the Hessian matrix estimated by finite-difference approximation.

One of the most widely used statistical models for tackling partially-bleached samples is the MAM (Galbraith et al., 1999; Galbraith and Roberts, 2012), which is based on well-established statistical principles. In this model, true logged characteristic D_e value is assumed to be drawn from a truncated Normal distribution, and the lower truncation point (i.e., γ) represents the mean of the logged characteristic dose of fully bleached grains. Two versions of the model are available, namely, the 3-parameter minimum age model (MAM₃) and its 4-parameter counterpart (MAM₄). The properties of the dataset (i.e., the degree of dispersion, the number of data points, the uncertainties associated with D_e values, etc.) may strongly affect the quality of the estimate in the application of the MAM₄ (Galbraith and Roberts, 2012). In contrast, the MAM₃ usually works well even for lowly dispersed D_e sets, while also being rather numerically stable compared to the MAM₄.

Galbraith et al. (1999) noted that some datasets, particularly those with small numbers of D_e values or less dispersed distributions, do not warrant fitting with the MAM₄, and that a more robust estimate of the characteristic dose may be obtained using the MAM₃. Consequently, the MAM₃ is conveniently regarded as the “default” model (Arnold et al., 2009) and has been widely applied in the literature to determine burial dose (e.g., Arnold et al., 2007; Schmidt et al., 2012; Kunz et al., 2013). According to Galbraith et al. (1999), the likelihood function of the MAM₃ can be described as follows: (2.2) p(Dj|θ)=MAM3(xj,yj|p,μ,σ)=p2πxj2e−(yj−μ)22xj2+2(1−p)[1−Φ(μ−μ0σ0)]2π(xj2+σ2)e(yj−μ)22(xj2+σ2) \begin{array}{*{35}{l}} p\left( {{D}_{j}}|\theta \right) & = & {{MAM}_{3}}\left( {{x}_{j}},{{y}_{j}}|p,\mu ,\sigma \right) \\ {} & = & \frac{p}{\sqrt{2\pi x_{j}^{2}}}{{e}^{-\frac{{{\left( {{y}_{j}}-\mu \right)}^{2}}}{2x_{j}^{2}}}}+\frac{2\left( 1-p \right)\left[ 1-\Phi \left( \frac{\mu -{{\mu }_{0}}}{{{\sigma }_{0}}} \right) \right]}{\sqrt{2\pi \left( x_{j}^{2}+{{\sigma }^{2}} \right)}}{{e}^{-\frac{{{\left( {{y}_{j}}-\mu \right)}^{2}}}{2\left( x_{j}^{2}+{{\sigma }^{2}} \right)}}} \end{array} where μ0=μσ2+yjxj21σ2+1xj2 {{\mu }_{0}}=\frac{\frac{\mu }{{{\sigma }^{2}}}+\frac{{{y}_{j}}}{x_{j}^{2}}}{\frac{1}{{{\sigma }^{2}}}+\frac{1}{x_{j}^{2}}} and σ0=11σ2+1xj2 {{\sigma }_{0}}=\frac{1}{\sqrt{\frac{1}{{{\sigma }^{2}}}+\frac{1}{x_{j}^{2}}}} .

Three quantities (p, μ, and σ) must be optimized in this model. p denotes the proportion of grains that are fully bleached before burial, μ denotes the logged characteristic dose, and σ is an adjustable parameter used to account for the dispersion in partially bleached grains (in unit per cent). Note that, conventionally, the logged characteristic dose in the MAM₃ is represented by γ, and in this study μ is alternatively used to make it correspond to that which is in the CAM. Φ(x) denotes the standard Normal cumulative distribution function. The optimal parameters are typically estimated by direct numerical maximization of the joint-likelihood that is determined using paired observations D={x_j, y_j | j=1,2,3,...,n). In contrast to the CAM, optimization of the MAM₃ is considered an ill-conditioned mathematical problem due to the strong nonlinearities in the calculation. Whether it is possible to converge it to global optimal estimates depends heavily on the choice of the initial parameters used. It is indispensable to try various initial parameter sets to obtain optimal estimates.

The maximum age model (MXAM) is an analogue of the MAM and is firstly applied by Olley et al. (2006) in their estimation of the maximum dose population from D_e distributions. Galbraith and Roberts (2012) enumerated two scenarios where the MXAM is applicable: (1) Samples composed of fully bleached grains that have subsequently been mixed with younger, intrusive grains, and (2) samples composed of fully bleached grains, some proportion of which have been bleached after deposition.

The characteristic dose of the MXAM₃ can be estimated according to the MAM₃ by using a “mirror image” of the original D_e distribution as an input (Olley et al., 2006; Galbraith and Roberts, 2012). Alternatively, by accounting for the fact that the lower truncation point of a truncated Normal distribution denotes the minimum characteristic dose in the MAM₃ while the upper truncation point of a truncated Normal distribution denotes the maximum characteristic dose in the MXAM₃, the likelihood for the MXAM₃ can be derived by a minor modification of Eq. 2.2 and can be described as follows: (2.3) p(Dj|θ)=MXAM3(xj,yj|p,μ,σ)=p2πxj2e−(yj−μ)22xj2+2(1−p)Φ(μ−μ0σ0)2π(xj2+σ2)e−(yj−μ)22(xj2+σ2) \begin{array}{*{35}{l}} p\left( {{D}_{j}}|\theta \right) & = & MXA{{M}_{3}}\left( {{x}_{j}},{{y}_{j}}|p,\mu ,\sigma \right) \\ {} & = & \frac{p}{\sqrt{2\pi x_{j}^{2}}}{{e}^{-\frac{{{\left( {{y}_{j}}-\mu \right)}^{2}}}{2x_{j}^{2}}}}+\frac{2\left( 1-p \right)\Phi \left( \frac{\mu -{{\mu }_{0}}}{{{\sigma }_{0}}} \right)}{\sqrt{2\pi \left( x_{j}^{2}+{{\sigma }^{2}} \right)}}{{e}^{-\frac{{{\left( {{y}_{j}}-\mu \right)}^{2}}}{2\left( x_{j}^{2}+{{\sigma }^{2}} \right)}}} \end{array} where μ₀ and σ₀ are calculated in the same manner as the corresponding values in the MAM₃.

Larger D_e estimate typically have larger absolute standard errors; namely, larger D_e estimates will vary more (Arnold and Roberts, 2009). Consequently, the frequency distribution of D_e estimates will tend to be positively skewed. The Lognormal distribution provides a simple model to describe variables wherein their standard deviations increase in proportion to their means. It is for this reason that the aforementioned statistical age models use log-transformed D_e datasets as inputs. Datasets in log-scale are obtained by calculating natural logarithms of D_e estimates and transforming their absolute standard errors into RSEs. This is because the RSE of a D_e estimate approximates the standard error of a D_e estimate on a natural logarithm scale (Galbraith and Roberts, 2012). As a result, these statistical age models are unsuitable for measured zero or negative D_e values that are consistent with zero or positive D_e values within a level of uncertainty. However, such types of D_e distribution are frequently encountered for either very young or modern samples. As a consequence, Arnold et al. (2009) used unlogged versions of the CAM and the MAM to analyze young/modern D_e datasets using a series of known-age fluvial samples and demonstrated that the unlogged versions of these models are capable of producing accurate burial dose estimates from their samples. To apply un-logged versions of these models, the user needs only use the original D_e values and their absolute standard errors as y_js and x_js, respectively, in Eqs. 2.1–2.3. In this case, the estimated μ denotes the characteristic dose (not the logged counterpart), while σ denotes a parameter in unit Gy (not in per cent).

The statistical age models described in Section 2 are in fact equivalent dose models. This is because they estimate the characteristic D_e instead of the burial age (e.g., Guérin et al., 2017). However, they can be easily transformed to obtain the posterior distribution of ages by taking the dose rate data into consideration during the simulation process: (3.1) μ=a×(d¯+ε) \mu =a\times \left( \bar{d}+\varepsilon \right) where a denotes the age, μ denotes the characteristic dose, d¯ \bar{d} denotes the mean annual dose rate, and ɛ denotes a variable that follows Gaussian distribution with mean 0 and standard deviation σ_d, namely, ɛ ∼ N(0, σ_d²). σ_d is calculated by aggregating many error sources due to measurements of uranium, thorium, potassium, water contents, etc. (Combès and Philippe, 2017). In this way, μ is represented mathematically using Eq. 3.1, and it is treated as an intermediate stochastic node during the simulation process. Consequently, the model can be re-parameterized using a instead of μ, namely, θ=[μ, σ] (in the CAM) and θ=[p, μ, σ] (in the MAM₃ or the MXAM₃) are re-parameterized as θ=[a, σ] and θ=[p, a, σ], respectively.

Each parameter is assigned a prior distribution that represents the prior knowledge of it before some evidence is considered. For p (i.e., the proportion of fully bleached grains), a Uniform prior that lies between 0.001 and 0.999 was used. For a (i.e., burial age) a typically employed non-informative prior is a Uniform distribution with lower and upper bounds set to equal a_min and a_max, which specifies a particular period of the burial history (e.g., Combès and Philippe, 2017). Parameter σ (i.e., a scale parameter inside the model) is an adjustable parameter that captures the unknown spread within the whole D_e distribution (for the CAM) or only the unknown spread in the D_e distribution for partially bleached grains (for the MAM₃). A significant amount of variation in OD exists between samples that either have different inherent luminescent properties, undergo different burial histories, and those that have been affected by various factors related to variation in microdosimetry and uncertainty from the measurement procedure. A summary of OD values from sedimentary samples (some are thought to have been fully bleached at the time of deposition) in published studies reveals that single-grain D_e distributions exhibit OD values of up to 30% (or higher), and the D_e distributions of multi-grain aliquots have a mean OD of 14% (Arnold and Roberts, 2009). OD values for several partially bleached young/modern fluvial samples calculated using an unlogged CAM range from 24% to 208% (Arnold et al., 2009). In this study, the prior of σ was represented using a Uniform distribution over the interval (0, 5).

Galbraith and Roberts (2012) reiterated that an estimate of additional uncertainty (σ_b) should be added in the quadrature to the RSE of each D_e value before running the MAM (as well as the MXAM) (Galbraith et al., 2005). σ_b represents the underlying spread in dose distribution typically present in well-bleached and undisturbed samples. An overestimate of σ_b will lead to an overestimate of the characteristic dose (and hence the burial age), and vice versa. Ideally, σ_b should be calculated from a well-bleached sample of the same mineral derived from the same source and, also ideally, of equivalent age (Galbraith and Roberts, 2012). In the absence of such information, σ_b values of 10% and 20% could be coarse approximations at multiple-grain and single-grain levels, respectively (e.g., Arnold and Roberts, 2009).

In this section, an age-depth model to combine D_e datasets from multiple samples with stratigraphic constraints was considered by using the statistical age models described in Section 2. Although Bayesian age-depth models are expected to improve the precision of dates, such precision may be misleading if potential systematic uncertainties have not been included by the individual samples (Zeeden et al., 2018). Compared to the application of statistical age models to individual samples using the standard approach, incorporating statistical age models into a Bayesian age-depth framework to refine the age-depth relationship between dates from multiple samples with order constraints requires that both random and systematic sources of errors are properly accounted for; otherwise, the precision of dates may increase decisively in unjustified ways (Combès and Philippe, 2017; Zeeden et al., 2018).

This study adopted the method described by Combès and Philippe (2017) to build the relationship between the characteristic dose μ_i, the annual dose rate d¯i {{\bar{d}}_{i}} , and the corresponding burial age a_i for the ith sample by taking into account both systematic and random sources of errors using the following Gaussian error model: (3.2) μi=ai×(d¯i+εdi+εdc) {{\mu }_{i}}={{a}_{i}}\times \left( {{{\bar{d}}}_{i}}+{{\varepsilon }_{di}}+{{\varepsilon }_{dc}} \right) where ɛ_di and ɛ_dc are the random measurement error for the ith sample and the systematic error shared by all samples from the sequence, respectively. The systematic error, which represents the error in the calibration of the measurement device (assuming that the measurements are made in the same luminescence laboratory), affects all estimates in the same way. The noise terms ɛ_di and ɛ_dc are assumed to follow Gaussian distributions, namely, εdi∼N(0,σdi2) {{\varepsilon }_{di}}\sim{\ }N\left( 0,\sigma _{di}^{2} \right) , and εdc∼N(0,σdc2) {{\varepsilon }_{dc}}\sim{\ }N\left( 0,\sigma _{dc}^{2} \right) , respectively. Therefore, the characteristic dose μ_i in the likelihoods is re-parameterized using a combination of burial age a and the dose rate data (i.e., the annual dose rate, associated random and systematic errors) according to Eq. 3.2 and can be treated as an intermediate stochastic node.

For parameters p_i and σ_i of the ith sample (i=1, 2, 3, ..., N), the priors were set equal to p and σ, respectively. For parameter a_i, a Uniform distribution for which the lower and upper bounds were set equal to a_mini and a_maxi was used (the same as Section 3 – MCMC sampling from statistical age models for individual samples). In addition, a constraint on vector [a₁, a₂, a₃, ..., a_N] such that a₁<a₂<a₃<...<a_N (i.e., wherein ages are sorted in ascending order) was imposed to allow stratigraphic information to be incorporated during the MCMC sampling process.

The objective being to obtain sampling distributions of parameters p(θ|D) based on the joint-likelihood p(D|θ) and priors p(θ) using MCMC sampling, a class of algorithms that utilize Markov chains to allow for the approximation of posterior distributions for parameters of interest by utilizing the prior distributions and likelihood as inputs (Annis et al., 2017). This study adopted the recently developed software package Stan (Gelman et al., 2015) for this purpose. Although there are dozens of built-in probability distributions that can be used to sample standard distributions (such as the Normal distribution, the Exponential distribution, etc.), users will occasionally require “non-standard” distributions (i.e., the joint-likelihoods for the CAM, MAM₃, and MXAM₃ used in this study) that may not yet be defined. Annis et al. (2017) gave a detailed tutorial on the implementation of MCMC sampling using user-defined distributions within the Stan programming language. Once the user-supplied function is defined, it can be called upon as a built-in distribution. This section provides a simple description of the numeric programs used to conduct MCMC sampling from statistical age models applying D_e datasets from both individual and multiple samples using the Stan software package. The relevant code and numeric scripts are freely downloadable from https://github.com/pengjunUCAS/mcmcSAM.

Individual samples

The input, likelihood, transformation, and distribution used to implement MCMC sampling for individual samples using the CAM and the MAM₃ are illustrated in Fig. 1. The MXAM₃ can be simulated in the same manner as the MAM₃ by changing the likelihood correspondingly. In Fig. 1, z_j and σ_zj denote the D_e and associated absolute error for the jth aliquot (grain). These statements can easily be modified to accommodate unlogged versions of the CAM, MAM₃, and MXAM₃. In Stan, variable definitions, distribution specifications, and program statements are placed within code blocks (Annis et al., 2017). Script “SAM0.stan” was employed to place all code blocks used to implement the sampling process. The R package “rstan” (an R interface for Stan) (Stan Development Team, 2018) was used to execute the script using a wrapper function mcmcSAM() available in the file “SAM.R”. Once the sampling process is terminated, the resultant MCMC convergence diagnostics are automatically stored in the file “mcmcSAM.pdf”, and the posterior samples are stored in the file ”mcmcSAM.csv”. Both files are available from the current R working directory. Please refer to the file “SAM.R” for more detail on arguments and default settings used in function mcmcSAM().

Fig. 1

Multiple samples with order constraints

Samples from the same sequence may substantially vary in D_e distributions due to differences in provenance, depositional environment, microdosimetry, etc. There is no generally applicable statistical age model that can be used to analyze all D_e distribution types. If we conceptualize three samples (S₁, S₂, and S₃) from the same stratigraphic sequence with their sample depths in ascending order, and we assume that these samples exhibit significantly different D_e distribution characteristics for which the most suitable models for these samples would be the MAM₃, CAM, and MXAM₃, respectively, then the statements used to conduct MCMC sampling on these samples applying a combination of different age models in log-scale can be summarized as shown in Fig. 2. z_ij and σ_zij denote the D_e and the associated absolute standard error for the jth aliquot (grain) from the ith sample. y_ij and x_ij denote the logged D_e and corresponding RSE for the jth aliquot (grain) from the ith sample. n₁, n₂, and n₃ denote the numbers of aliquots (grains) for the three samples. p(D|θ) is calculated as the product of joint likelihoods for samples S₁, S₂, and S₃ that are fitted using the MAM₃, CAM, and MXAM₃, respectively. Note that in this study only three samples are used for illustrative purposes, while in practice many samples fitted using various combinations of the three different age models can be simultaneously analyzed in a manner similar to that demonstrated in Fig. 2. This allows the programs to have a degree of flexibility when D_e distributions from the same sequence need to be treated differently. Two scripts (i.e., “SAMseq0.stan” and “SAMseq1.stan”) were used to conduct the sampling process without and with order constraints, respectively. A wrapper function mcmcSAMseq()that is available in file “SAMseq.R” was used to call these scripts. The resultant MCMC convergence diagnostics are automatically stored in file “mcmcSAMseq.pdf”, and the posterior samples are stored in file “mcmcSAMseq.csv”.

Fig. 2

First, the validity of the MCMC sampling protocol (as demonstrated in Fig. 1) used to analyze D_e sets from individual samples was checked using measured multi-grain D_e distributions (Fig. 3) of four aeolian samples (i.e., GL2-1 to GL2-4) from Gulang County at the southern margin of the Tengger Desert, China (Peng et al., 2016b). The D_e distributions were obtained using the single aliquot regenerative-dose (SAR) protocol (Galbraith et al., 1999; Murray and Wintle, 2000). When applying the MAM₃ and MXAM₃, an additional uncertainty (i.e., σ_b) of 10% was added in the quadrature to the RSEs of the measured D_e values to account for unrecognized measurement errors. In this study (unless stated otherwise), three parallel Markov chains (the number of iterations for each chain was set equal to 10000) were simulated during the MCMC sampling process for the individual samples, and four parallel Markov chains (the number of iterations for each chain was set equal to 6000) were simulated during the MCMC sampling process for the multiple samples. The lower and upper bounds for the priors of ages were set equal to a_min(a_mini)=1 ka and a_max(a_maxi)=100 ka, respectively. The characteristic dose (i.e., μ) and age (i.e., a) of the CAM, MAM₃, and MXAM₃ obtained using MLE and MCMC were then compared and the results summarized in Table 1. The estimates obtained between the two different methods were consistent within errors, suggesting that estimates derived from MCMC can be used as an informative comparison for those estimated by MLE.

Fig. 3

These four samples were also analyzed simultaneously using the MCMC sampling protocol (as demonstrated in Fig. 2) for multiple samples. D_e values of these samples fell within reasonably narrow ranges and displayed homogeneous distributions (Fig. 3). Accordingly, the CAM was applied to determine the burial dose. Samples GL2-1, GL2-2, GL2-3, and GL2-4 were collected from the same section in ascending order depths. Consequently, their burial ages were also expected to be in ascending order. The MCMC sampling protocol was applied to these samples with order constraints (results shown in Table 2). The systematic error representing the error in the calibration of the dose rate measurement device shared by all samples was assumed to be σ_dc=0.1. The posterior distributions of burial age are shown in Fig. 4. Results obtained without order constraints were also provided for comparison. These results suggested that burial ages would be internally consistent with the stratigraphy if order constraints were imposed during the MCMC sampling process. Moreover, the posterior standard deviations of burial ages would be significantly reduced when using stratigraphic constraints, implying that the precision of burial ages improved. Finally, a comparison between Table 1 and Table 2 for CAM results demonstrated that burial dose and ages (and their posterior standard deviations) obtained by MCMC sampling using both individual samples and multiple samples simultaneously without order constraints were indistinguishable from each other.

Fig. 4

This section shows simulation results of a series of D_e distributions with known burial ages to further verify the performance of the MCMC sampling protocols. The first three samples (#1–#3) were assumed to be partially-bleached, the next three samples (#4–#6) fully-bleached, and the last three samples (#7–#9) contain fully-bleached grains that have been subsequently mixed with younger, intrusive grains. These samples were assumed to be collected from the same sedimentary sequence but at different depths. Their “true” burial doses ranged from 18 to 42 Gy with a step of 3 Gy, and their “true” burial ages ranged from 7 to 11 ka with a step of 0.5 ka. The “true” average dose rate can be obtained by dividing the “true” burial dose by the “true” burial age of a sample. The random error arising from dose rate measurement was assumed to be 7%, in order to simulate a realistic annual dose rate. Namely, the “true” dose rates are noised using a Normal distribution with a mean equivalent to the “true” average dose rate and a relative standard deviation (RSD) equal to 7%. The generated random dose rates were used as the measured dose rates when determining burial ages using MLE and MCMC. The systematic error representing the error in the calibration of the dose rate measurement device shared by all samples was assumed to be zero (i.e., σ_dc=0).

Lognormal distributions were used to model the natural dispersion in dose distributions arising from dose rate variation with an RSD of 5% (i.e., σ_d=0.05). Measured D_e distributions were simulated according to the method described by Li et al. (2017), which involves simulating OSL intensities of the natural dose (L_n), a series of regenerative doses (L_x), and their corresponding test doses (T_n and T_x). Measured OSL sensitivity data taken from an empirical sample was used as the basis of the simulation. The OSL sensitivity (e.g., counts per unit of dose) of the different grains were randomly generated from a Gamma distribution fitted to the single-grain T_n datum from a sample taken from Lake Mungo, Australia (Fig. 5). Noise was added to each of the simulated OSL signal intensities based on the uncertainty arising from counting statistics and instrumental irreproducibility. It was assumed that both photon counts and dark counts were over-dispersed compared to a Poisson distribution (e.g., Li, 2007; Adamiec et al., 2012) and follow Negative Binomial distributions (e.g., Bluszcz et al., 2015). The variance to mean ratios for the photon counts and dark counts were set equal to 1.38 and 1.85, respectively. The dark count was set equal to 10 cts/0.2s. The instrumental reproducibility was set as 2%. For the natural signals (L_n), an extra noise from “intrinsic OD” (e.g., σ_m=10%), which represents unrecognized measurement errors (e.g., OD observed from a dose recovery experiment) (e.g., Thomsen et al., 2005; Jacobs et al., 2006; Guérin et al., 2017), was also added. D_e estimation was conducted using the R package “numOSL” (Peng et al., 2013; Peng and Li, 2017). For each sample, 100 D_e values (and associated errors) were simulated. The nine simulated D_e distributions are shown in Fig. 6.

Fig. 5

Fig. 6

The simulated D_e sets were analyzed individually using MLE, analyzed simultaneously using MCMC without order constraints, and analyzed simultaneously using MCMC with order constraints. The first, next, and last three samples were analyzed using the MAM₃, CAM, and MXAM₃, respectively. To facilitate the model in the analysis of the MAM₃ and the MXAM₃, a σ_b should be determined beforehand. The calculated OD values for samples #4–#6 can be regarded as appropriate σ_b values, given that these samples are “fully-bleached” simulated samples whose errors, arising from dose rate variation (i.e., σ_d) and unrecognized measurement uncertainty (i.e., σ_m), were the same as samples #1–#3 and samples #7–#9. The estimated OD values were 10.09±0.31%, 10.48±0.33%, and 11.89±0.38%, respectively, for samples #4–#6. These values are broadly consistent with the expected value of σd2+σm2=(5%)2+(10%)2=11.18% \sqrt{\sigma _{d}^{2}+\sigma _{m}^{2}}=\sqrt{{{\left( 5\% \right)}^{2}}+{{\left( 10\% \right)}^{2}}}=11.18\% , confirming the validity of the simulation. Accordingly, a σ_b value of 11% was added to the RSEs of the D_e values when analyzing the MAM₃ and the MXAM₃. Analyzed results for each sample are summarized in Table 3. The posterior distributions of burial ages are shown in Fig. 7 for results obtained by MCMC sampling with and without order constrains. Similar to Section 4 (Measured D_e sets), the burial doses and ages (and associated errors) obtained by MLE using individual samples and by MCMC using multiple samples simultaneously without order constraints were broadly consistent between one another. When order constraints were imposed during the MCMC sampling process, the precision of burial age improved for all samples, and the accuracy of burial age also increased for most samples (except for samples #4 and #9) (Fig. 7). The relative error (RE) of burial age of sample #4 increased slightly from 2.09% to 2.20%, and that of sample #9 increased substantially from 0.35% to 4.4%, if order constraints were imposed.

Fig. 7

Although a properly constructed MCMC protocol enables one to draw samples successively from a convergent Markov chain, deciding the point to terminate sampling is sometimes difficult as well as the point to confidently conclude when the algorithm has converged to the desired stationary distribution. One simple and direct way to inspect the convergence of a chain is to check the degree of mixing of the chain by using a trace plot (Fig. 8A). Values that get “stuck” within certain space intervals of variables indicate poor mixing. In contrast, a chain that shows good mixing properties will move freely along the feasible space of variables. Another simple method for assessing the convergence of a chain is to monitor the autocorrelations of generated variables (Fig. 8C). Variables generated using MCMC suffer from autocorrelations to some extent. To decrease the autocorrelations of a variable, simulated samplers are routinely trimmed using a “thinning” protocol after the “burn-in” (“warm-up”) procedure. If the autocorrelations of variables remain extremely high after applying a large “thinning” value, it may imply that the speed of the convergence is slow (i.e., poor mixing) and more samples are therefore needed so that a meaningful inference can be drawn. In this study (unless stated otherwise), the number of “warm-up” iterations per chain was set equal to 2000, and the number of “thinning” iterations per chain was set equal to 1.

Fig. 8

The Gelman-Rubin convergence diagnostic (Gelman and Rubin, 1992) is one of the most widely used methods for monitoring the convergence of MCMC outputs. Several parallel chains are simulated using various initial states when applying this method, and a shrink factor is used as a measurement of the difference between within-chain and between-chain variance (which is similar to the analysis of variance). A shrink factor value equal to or less than 1 indicates adequate convergence, while a factor value significantly above 1 indicates a lack of convergence. Results from Gelman-Rubin convergence diagnostics for posterior samples of CAM burial ages from sample GL2-1 is shown in Fig. 8D. It can be clearly observed that the simulations stabilized and were very close to 1 after approximately 3000 iterations. A summary of Gelman-Rubin convergence diagnostics for sequences of measured and simulated samples obtained using the MCMC sampling protocol of Fig. 2 with order constraints is presented in Table 4.