In optically stimulated luminescence (OSL) dating, statistical age models for equivalent dose (D_{e}) distributions are routinely estimated using the maximum likelihood estimation (MLE) method. In this study, a Markov chain Monte Carlo (MCMC) method was used to analyze statistical age models, including the central age model (CAM), the minimum age model (MAM), the maximum age model (MXAM), _{e} distributions from individual sedimentary samples and subsequently extended to simultaneously extract age estimates from multiple samples with stratigraphic constraints. The MCMC method allows for the use of Bayesian inference to refine chronological sequences from multiple samples, including both fully and partially bleached OSL dates. This study designed easily implemented open-source numeric programs to perform MCMC sampling. Measured and simulated D_{e} distributions are used to validate the reliability of dose (age) estimates obtained by this method. Findings from this study demonstrate that estimates obtained by the MCMC method can be used to informatively compare results obtained by the MLE method. The application of statistical age models to multiple OSL dates with stratigraphic orders using the MCMC method may significantly improve both the precision and accuracy of burial ages.

#### Keywords

- OSL dating
- statistical age models
- equivalent dose
- MCMC sampling

In optically stimulated luminescence (OSL) dating, the assessment of certain parameters of interest (_{e}) is essential in estimating sample age. Routinely, a large set of aliquots (grains) are analyzed from which a single D_{e} (_{e} determined using a statistical model) is ultimately obtained along with its standard error estimate. Galbraith and Roberts (2012) reviewed the most commonly used statistical age models in determining the characteristic D_{e}, namely, from the simple central age model (CAM) (Galbraith

Maximum likelihood estimation (MLE), which provides point estimates of parameters and associated standard errors, is routinely adopted to estimate parameters used in these statistical age models by maximizing the logged-likelihood function (Galbraith _{e} estimated using MLE by the average environmental dose rate of the sample, and the associated standard error is estimated using propagation of error formulas (referred to as the standard approach). However, the standard approach only provides a coarse characterization of the statistics associated with burial age. Moreover, it does not permit the inclusion of additional information, such as order constraints between the ages of a given stratigraphic sequence to improve age estimates (Combès and Philippe, 2017).

In contrast, the Markov chain Monte Carlo (MCMC) methods construct a Markov chain while converging to the desired sampling distribution by starting from random initials and seeking to characterize posterior distributions for parameters of interest (Gelman

Such MCMC methods have increasingly been employed to analyze OSL data for Bayesian inferences (e.g., Millard, 2004; Sivia _{e} sets. Cunningham _{e} distributions while using an MCMC method to obtain the posterior of parameter estimates. Combès _{e} for fully bleached aliquots (grains). A Metropolis-within-Gibbs sampler is used to obtain the posterior of the characteristic dose. The model has been further expanded to analyse poorly bleached samples using a Gaussian mixture model (e.g., Christophe

The first objective of this study is to estimate statistical model parameters (such as the CAM, the MAM,

Bayesian models utilize stratigraphic relationships between samples to increase the precision of dates and to provide a sophisticated protocol to combine multiple age estimates into a meaningful chronological framework (Ramsey, 1995, 2008). By incorporating the stratigraphic principle that shallower sediments in the stratum sequence are younger than deeper sediments into age-depth models, Bayesian inference can reduce the uncertainty of dates wherein age probability distributions overlap. Age-depth modelling using Bayesian inference was initially developed to analyze radiocarbon dating data (Ramsey, 1995; Bayliss and Ramsey, 2004). Over the past decade, the application of Bayesian inference to model depth-age relationships using luminescence data have also gained widespread usage (e.g., Rhodes

During the routine application of statistical age models, such as the CAM and MAM, characteristic D_{e} values of individual samples are estimated separately (e.g., Schmidt _{e} for a sequence of young fluvial samples using an unlogged version of the MAM. The replicate characteristic D_{e} values obtained by bootstrapping are used to construct a probability density function for each sample in order to apply Bayesian inference to improve OSL chronologies. However, as reiterated by Cunningham and Wallinga (2012), the boot-strap likelihood protocol is just an analogue of the likelihood function and does not produce a true likelihood.

The second objective of this study is to apply statistical age models to refine age-depth relationships using D_{e} datasets from multiple luminescence samples employing an alternative Bayesian inference model. In contrast to the method by Cunningham and Wallinga (2012), in this study characteristic D_{e} values are sampled directly from the probability density function that describes the statistical model (see the following section for more detail) through MCMC sampling. The joint probability density function for multiple OSL ages is constructed by utilizing a combination of statistical models and measured D_{e} and dose rate datasets. Both the CAM and MAM are used to flexibly analyze fully and partially bleached OSL dates, respectively, using MCMC in order to combine different types of statistical age models in a mathematically rigorous, geologically reasonable, and philosophically satisfying manner. The performance of the model is then tested on measured and simulated D_{e} distributions.

Galbraith _{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 _{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:
_{j}_{j}_{e} value and its relative standard error (RSE), respectively. p(_{j}_{j}_{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 _{j}_{j}

One of the most widely used statistical models for tackling partially-bleached samples is the MAM (Galbraith _{e} value is assumed to be drawn from a truncated Normal distribution, and the lower truncation point (_{3}) and its 4-parameter counterpart (MAM_{4}). The properties of the dataset (_{e} values, _{4} (Galbraith and Roberts, 2012). In contrast, the MAM_{3} usually works well even for lowly dispersed D_{e} sets, while also being rather numerically stable compared to the MAM_{4}.

Galbraith _{e} values or less dispersed distributions, do not warrant fitting with the MAM_{4}, and that a more robust estimate of the characteristic dose may be obtained using the MAM_{3}. Consequently, the MAM_{3} is conveniently regarded as the “default” model (Arnold _{3} can be described as follows:

Three quantities (_{3} is represented by _{j}_{j}_{3} 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 _{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_{3} can be estimated according to the MAM_{3} by using a “mirror image” of the original D_{e} distribution as an input (Olley _{3} while the upper truncation point of a truncated Normal distribution denotes the maximum characteristic dose in the MXAM_{3}, the likelihood for the MXAM_{3} can be derived by a minor modification of _{0} and _{0} are calculated in the same manner as the corresponding values in the MAM_{3}.

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 _{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_{j}s and x_{j}s, respectively, in

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 _{d}_{d}^{2}). _{d}_{3} or the MXAM_{3}) are re-parameterized as

Each parameter is assigned a prior distribution that represents the prior knowledge of it before some evidence is considered. For _{min}_{max}_{e} distribution (for the CAM) or only the unknown spread in the D_{e} distribution for partially bleached grains (for the MAM_{3}). 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

Galbraith and Roberts (2012) reiterated that an estimate of additional uncertainty (_{b}_{e} value before running the MAM (as well as the MXAM) (Galbraith _{b}_{b}_{b}_{b}

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

This study adopted the method described by Combès and Philippe (2017) to build the relationship between the characteristic dose _{i}_{i}_{di}_{dc}_{di} and _{dc} are assumed to follow Gaussian distributions, namely,
_{i}

For parameters _{i} and _{i}_{i}, a Uniform distribution for which the lower and upper bounds were set equal to _{mini}_{maxi}_{1}, _{2}, _{3}, ..., _{N}] such that _{1}<_{2}<_{3}<...<_{N} (

The objective being to obtain sampling distributions of parameters p(θ|D) based on the joint-likelihood p(_{3}, and MXAM_{3} used in this study) that may not yet be defined. Annis _{e} datasets from both individual and multiple samples using the Stan software package. The relevant code and numeric scripts are freely downloadable from

The input, likelihood, transformation, and distribution used to implement MCMC sampling for individual samples using the CAM and the MAM_{3} are illustrated in _{3} can be simulated in the same manner as the MAM_{3} by changing the likelihood correspondingly. In _{j}_{zj}_{e} and associated absolute error for the _{3}, and MXAM_{3}. In Stan, variable definitions, distribution specifications, and program statements are placed within code blocks (Annis

Samples from the same sequence may substantially vary in D_{e} distributions due to differences in provenance, depositional environment, microdosimetry, _{e} distribution types. If we conceptualize three samples (S_{1}, S_{2}, and S_{3}) 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_{3}, CAM, and MXAM_{3}, 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 _{ij}_{zij}_{e} and the associated absolute standard error for the _{ij}_{ij}_{e} and corresponding RSE for the _{1}, _{2}, and _{3} denote the numbers of aliquots (grains) for the three samples. _{1}, S_{2}, and S_{3} that are fitted using the MAM_{3}, CAM, and MXAM_{3}, 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 _{e} distributions from the same sequence need to be treated differently. Two scripts (

_{e}sets

First, the validity of the MCMC sampling protocol (as demonstrated in _{e} sets from individual samples was checked using measured multi-grain D_{e} distributions (_{e} distributions were obtained using the single aliquot regenerative-dose (SAR) protocol (Galbraith _{3} and MXAM_{3}, an additional uncertainty (_{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 _{min}_{mini}_{max}_{maxi}_{3}, and MXAM_{3} obtained using MLE and MCMC were then compared and the results summarized in

Comparisons of burial dose and age estimated from MLE and MCMC for the various statistical age models using D_{e} sets from four measured aeolian samples. Quantities estimated using the MCMC sampling protocol shown in Fig. 1 are marked in bold. Quantities estimated using the MLE are inside parentheses.

_{3} |
_{3} |
||||||
---|---|---|---|---|---|---|---|

GL2-1 | 3.26 ± 0.25 | ||||||

GL2-2 | 2.84 ± 0.21 | ||||||

GL2-3 | 3.23 ± 0.23 | ||||||

GL2-4 | 3.40 ± 0.24 |

These four samples were also analyzed simultaneously using the MCMC sampling protocol (as demonstrated in _{e} values of these samples fell within reasonably narrow ranges and displayed homogeneous distributions (_{dc}=0.1. The posterior distributions of burial age are shown in

A summary of burial dose and age estimated from the CAM for D_{e} sets from four measured aeolian samples using the MCMC sampling protocol shown in Fig. 2 without and with order constraints. The systematic error for the dose rate measurement shared by all samples was set to σ_{dc}=0.1.

GL2-1 | 3.26 ± 0.25 | 34.55 ± 0.91 | 10.75 ± 0.96 | 34.23 ± 0.86 | 9.72 ± 0.52 |

GL2-2 | 2.84 ± 0.21 | 32.71 ± 1.00 | 11.68 ± 1.06 | 32.19 ± 0.92 | 10.28 ± 0.49 |

GL2-3 | 3.23 ± 0.23 | 30.58 ± 1.04 | 9.59 ± 0.85 | 31.03 ± 1.00 | 10.59 ± 0.52 |

GL2-4 | 3.40 ± 0.24 | 32.24 ± 1.29 | 9.60 ± 0.85 | 33.13 ± 1.26 | 11.13 ± 0.67 |

_{e}sets

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 (_{dc}=0).

Lognormal distributions were used to model the natural dispersion in dose distributions arising from dose rate variation with an RSD of 5% (_{d}=0.05). Measured D_{e} distributions were simulated according to the method described by Li _{n}), a series of regenerative doses (_{x}), and their corresponding test doses (_{n} and _{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 _{n} datum from a sample taken from Lake Mungo, Australia (_{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 _{e} estimation was conducted using the R package “numOSL” (Peng _{e} values (and associated errors) were simulated. The nine simulated D_{e} distributions are shown in

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_{3}, CAM, and MXAM_{3}, respectively. To facilitate the model in the analysis of the MAM_{3} and the MXAM_{3}, 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 (_{d}) and unrecognized measurement uncertainty (_{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
_{b} value of 11% was added to the RSEs of the D_{e} values when analyzing the MAM_{3} and the MXAM_{3}. Analyzed results for each sample are summarized in _{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) (

A summary of burial dose and age estimated from various statistical age models for nine simulated D_{e} sets taken from the same sedimentary sequence. Results obtained using MLE and results obtained using the MCMC sampling protocol as shown in Fig. 2 without and with order constraints are presented. The systematic error on the dose rate measurement shared by all simulated samples was assumed to be σ_{dc}=0.

#1 | MAM_{3} |
7 | 18 | 18.35 ± 1.17 | 7.95 ± 0.75 | 17.85 ± 1.17 | 7.82 ± 0.76 | 16.58 ± 1.02 | 6.80 ± 0.44 |

#2 | 7.5 | 21 | 21.16 ± 0.63 | 6.95 ± 0.53 | 20.96 ± 0.79 | 6.96 ± 0.57 | 21.17 ± 0.65 | 7.23 ± 0.40 | |

#3 | 8 | 24 | 24.47 ± 1.22 | 8.41 ± 0.72 | 23.64 ± 1.57 | 8.22 ± 0.81 | 23.25 ± 1.43 | 7.84 ± 0.42 | |

#4 | CAM | 8.5 | 27 | 26.91 ± 0.33 | 8.24 ± 0.59 | 26.92 ± 0.34 | 8.32 ± 0.61 | 26.93 ± 0.33 | 8.31 ± 0.39 |

#5 | 9 | 30 | 30.40 ± 0.38 | 8.50 ± 0.60 | 30.41 ± 0.39 | 8.59 ± 0.63 | 30.44 ± 0.38 | 8.83 ± 0.42 | |

#6 | 9.5 | 33 | 32.79 ± 0.46 | 10.59 ± 0.76 | 32.81 ± 0.47 | 10.71 ± 0.79 | 32.66 ± 0.45 | 9.62 ± 0.42 | |

#7 | MXAM_{3} |
10 | 36 | 35.77 ± 1.61 | 9.13 ± 0.76 | 36.43 ± 1.86 | 9.39 ± 0.84 | 37.18 ± 1.81 | 9.95 ± 0.45 |

#8 | 10.5 | 39 | 39.26 ± 1.37 | 9.32 ± 0.73 | 39.70 ± 1.69 | 9.53 ± 0.79 | 40.62 ± 1.87 | 10.39 ± 0.54 | |

#9 | 11 | 42 | 43.65 ± 2.55 | 10.63 ± 0.97 | 44.61 ± 2.52 | 10.96 ± 1.00 | 45.43 ± 2.53 | 11.48 ± 0.85 |

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 (

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

A summary of Gelman-Rubin convergence diagnostics for measured and simulated samples obtained using the MCMC sampling protocol, as shown in Fig. 2 with order constraints. n_eff is the effective sample size, while Rhat (i.e., the shrink factor) is a statistic measure of the ratio of the average variance of samples within each chain to the variance of the pooled samples across chains. If all chains are at equilibrium, the Rhat will be 1. If these chains have not converged to a common distribution, the Rhat statistic will be greater than 1.

GL2-1 | 16000 | 0.9998 | 5541.74 | 1.0002 | #1 | 13865.37 | 1.0001 | 9987.93 | 1.0000 |

GL2-2 | 16000 | 0.9999 | 8220.18 | 1.0003 | #2 | 16000 | 1.0001 | 11473.55 | 1.0000 |

GL2-3 | 16000 | 0.9997 | 7595.60 | 1.0002 | #3 | 9405.95 | 1.0004 | 12466.84 | 1.0003 |

GL2-4 | 16000 | 0.9998 | 8864.01 | 1.0003 | #4 | 16000 | 1.0003 | 13161.07 | 0.9999 |

#5 | 16000 | 0.9999 | 12830.44 | 0.9999 | |||||

#6 | 16000 | 0.9999 | 14048.98 | 0.9998 | |||||

#7 | 16000 | 0.9999 | 13073.69 | 0.9998 | |||||

#8 | 7780.82 | 1.0003 | 12665.56 | 1.0001 | |||||

#9 | 16000 | 0.9999 | 16000 | 0.9999 |

This study adopted an MCMC method to obtain sampling distributions on parameters of interest in statistical age models, including the MAM_{3}, CAM, and MXAM_{3}, using D_{e} distributions from individual samples and multiple samples with order constraints. It demonstrated that the estimates obtained by MLE and MCMC from measured samples were consistent within errors for various age models (

For all MCMC simulations, initial states did not need to be specified and were generated automatically using random number generators, suggesting that the MCMC sampling method employed in this study was not sensitive to the choice of the initial state, and that the algorithm itself was free from the influence of local modes and was able to converge to the desired distribution with a finite number of iterations (_{3}, and the MXAM_{3} (Peng

When applying the MCMC sampling method to multiple samples with stratigraphic constraints using statistical age models, samples are not analyzed on a standalone basis; rather, they are combined simultaneously in a single model implementation so that between-sample age-depth relationships can also be included. This shows that the burial age of measured samples was internally consistent with the stratigraphy, and the posterior standard deviations of burial ages significantly decreased (_{e} sets with known burial ages demonstrate that the MCMC sampling method can not only increase the precision but also may improve the accuracy of age estimates (_{e} sets from multiple samples (e.g., Cunningham and Wallinga, 2012; Combès and Philippe, 2017). However, it should be noted that incorporating age-depth relationship between samples into the Bayesian model may also have the risk of decreasing the accuracy of age estimates for certain samples within the sequence, as observed in

The sequence of D_{e} sets simulated in Section 4 (Simulated D_{e} sets) was determined to be free of systematic error arising from dose rate measurements. However, the systematic error term that induces dependences between observations and the associated age-depth model should always be taken into account when analyzing measured samples (e.g., Rhodes _{dc}=0.1 for measured samples in this study. The quality of estimates for burial ages may potentially be improved by introducing external (independent) dates, allowing for the correction of systematic components of the error (e.g., Rhodes

This study employed an MCMC sampling method for the statistical analysis of D_{e} datasets, using statistical age models, including the CAM, MAM_{3}, and MXAM_{3}. The method was tested using measured, and simulated D_{e} sets, and consistent results in agreement with expected values were obtained. The development of numerical programs provides a statistical toolbox for the calculation of OSL ages using MCMC sampling for both individual samples and multiple samples with (or without) order constraints by taking into account systematic and individual errors. The MCMC method employed in this study can be used in addition to other Bayesian models in tackling luminescence-related chronological data and can be flexibly adapted to compute OSL ages of both well- and poorly-bleached sedimentary samples.

#### A summary of burial dose and age estimated from various statistical age models for nine simulated De sets taken from the same sedimentary sequence. Results obtained using MLE and results obtained using the MCMC sampling protocol as shown in Fig. 2 without and with order constraints are presented. The systematic error on the dose rate measurement shared by all simulated samples was assumed to be σdc=0.

#1 | MAM_{3} |
7 | 18 | 18.35 ± 1.17 | 7.95 ± 0.75 | 17.85 ± 1.17 | 7.82 ± 0.76 | 16.58 ± 1.02 | 6.80 ± 0.44 |

#2 | 7.5 | 21 | 21.16 ± 0.63 | 6.95 ± 0.53 | 20.96 ± 0.79 | 6.96 ± 0.57 | 21.17 ± 0.65 | 7.23 ± 0.40 | |

#3 | 8 | 24 | 24.47 ± 1.22 | 8.41 ± 0.72 | 23.64 ± 1.57 | 8.22 ± 0.81 | 23.25 ± 1.43 | 7.84 ± 0.42 | |

#4 | CAM | 8.5 | 27 | 26.91 ± 0.33 | 8.24 ± 0.59 | 26.92 ± 0.34 | 8.32 ± 0.61 | 26.93 ± 0.33 | 8.31 ± 0.39 |

#5 | 9 | 30 | 30.40 ± 0.38 | 8.50 ± 0.60 | 30.41 ± 0.39 | 8.59 ± 0.63 | 30.44 ± 0.38 | 8.83 ± 0.42 | |

#6 | 9.5 | 33 | 32.79 ± 0.46 | 10.59 ± 0.76 | 32.81 ± 0.47 | 10.71 ± 0.79 | 32.66 ± 0.45 | 9.62 ± 0.42 | |

#7 | MXAM_{3} |
10 | 36 | 35.77 ± 1.61 | 9.13 ± 0.76 | 36.43 ± 1.86 | 9.39 ± 0.84 | 37.18 ± 1.81 | 9.95 ± 0.45 |

#8 | 10.5 | 39 | 39.26 ± 1.37 | 9.32 ± 0.73 | 39.70 ± 1.69 | 9.53 ± 0.79 | 40.62 ± 1.87 | 10.39 ± 0.54 | |

#9 | 11 | 42 | 43.65 ± 2.55 | 10.63 ± 0.97 | 44.61 ± 2.52 | 10.96 ± 1.00 | 45.43 ± 2.53 | 11.48 ± 0.85 |

#### A summary of Gelman-Rubin convergence diagnostics for measured and simulated samples obtained using the MCMC sampling protocol, as shown in Fig. 2 with order constraints. n.eff is the effective sample size, while Rhat (i.e., the shrink factor) is a statistic measure of the ratio of the average variance of samples within each chain to the variance of the pooled samples across chains. If all chains are at equilibrium, the Rhat will be 1. If these chains have not converged to a common distribution, the Rhat statistic will be greater than 1.

GL2-1 | 16000 | 0.9998 | 5541.74 | 1.0002 | #1 | 13865.37 | 1.0001 | 9987.93 | 1.0000 |

GL2-2 | 16000 | 0.9999 | 8220.18 | 1.0003 | #2 | 16000 | 1.0001 | 11473.55 | 1.0000 |

GL2-3 | 16000 | 0.9997 | 7595.60 | 1.0002 | #3 | 9405.95 | 1.0004 | 12466.84 | 1.0003 |

GL2-4 | 16000 | 0.9998 | 8864.01 | 1.0003 | #4 | 16000 | 1.0003 | 13161.07 | 0.9999 |

#5 | 16000 | 0.9999 | 12830.44 | 0.9999 | |||||

#6 | 16000 | 0.9999 | 14048.98 | 0.9998 | |||||

#7 | 16000 | 0.9999 | 13073.69 | 0.9998 | |||||

#8 | 7780.82 | 1.0003 | 12665.56 | 1.0001 | |||||

#9 | 16000 | 0.9999 | 16000 | 0.9999 |

#### A summary of burial dose and age estimated from the CAM for De sets from four measured aeolian samples using the MCMC sampling protocol shown in Fig. 2 without and with order constraints. The systematic error for the dose rate measurement shared by all samples was set to σdc=0.1.

GL2-1 | 3.26 ± 0.25 | 34.55 ± 0.91 | 10.75 ± 0.96 | 34.23 ± 0.86 | 9.72 ± 0.52 |

GL2-2 | 2.84 ± 0.21 | 32.71 ± 1.00 | 11.68 ± 1.06 | 32.19 ± 0.92 | 10.28 ± 0.49 |

GL2-3 | 3.23 ± 0.23 | 30.58 ± 1.04 | 9.59 ± 0.85 | 31.03 ± 1.00 | 10.59 ± 0.52 |

GL2-4 | 3.40 ± 0.24 | 32.24 ± 1.29 | 9.60 ± 0.85 | 33.13 ± 1.26 | 11.13 ± 0.67 |

#### Comparisons of burial dose and age estimated from MLE and MCMC for the various statistical age models using De sets from four measured aeolian samples. Quantities estimated using the MCMC sampling protocol shown in Fig. 1 are marked in bold. Quantities estimated using the MLE are inside parentheses.

_{3} |
_{3} |
||||||
---|---|---|---|---|---|---|---|

GL2-1 | 3.26 ± 0.25 | ||||||

GL2-2 | 2.84 ± 0.21 | ||||||

GL2-3 | 3.23 ± 0.23 | ||||||

GL2-4 | 3.40 ± 0.24 |

_{e}) distributions: Implications for OSL dating of sediment mixtures. _{e}) distributions: Implications for OSL dating of sediment mixtures

^{2}). ^{2})

_{e} distributions and calculation of D_{e}. _{e} distributions and calculation of D_{e}