Philippe Lanos

Anne Philippe

HIERARCHICAL BAYESIAN MODELING FOR COMBINING DATES IN ARCHAEOLOGICAL CONTEXT

Keywords: Bayesian statistics, ChronoModel software, individual errors, MCMC methods, outlier penalization

A Bayesian approach is proposed for combining dates from different dating methods used in archeology. This modeling provides an automatic way to penalize outlying data. Examples are provided from different archeological contexts and involving radiocarbon, luminescence and archaeomagnetic results. This new combination procedure is also applied to the wiggle-matching process in dendrochronological dating. Calculations are based on MCMC numerical techniques and can be performed using the cross-platform ChronoModel application which is free, open source software (FOSS).

Introduction

Many dating methods are available to learn about the age of archeological artifacts. The choice of the method generally depends on the nature of the object and its age (see [START_REF] Aitken | Science Based Dating in Archaeology[END_REF] for a detailed review). Radiocarbon, luminescence, and archaeomagnetic dating are probably the most commonly used methods.

Radiocarbon dating (or 14C dating) is used to evaluate the ages of biological specimens (e.g. charcoals, wood artifacts, bone remains, etc.). It is possible to date objects as old as around 50 000 years. Many statistical papers are devoted to this method and software applications are available to analyze and model the radiocarbon data (see [START_REF] Buck | BCal : an on-line Bayesian radiocarbon calibration tool[END_REF][START_REF] Bronk Ramsey | Recent and planned developments of the program OxCal[END_REF]Bronk Ramsey, 2009a).

Luminescence dating (TL/OSL) is based on the radiation absorbed and stored in the crystal lattice of minerals such as quartz. The calculated age is the elapsed time since the last exposure to sunlight or intense heat.

Archaeomagnetic dating (AM) relies on the past variations of the Earth's magnetic field. The ferromagnetic minerals, heated to high temperatures, acquire a remnant magnetization with a direction parallel to, and an intensity proportional to, the local Earth's magnetic field at the time of cooling. The method is implemented for some geographical area where variations of the Earth's magnetic field are known. It allows the last heating of baked clays to be dated.

One important question in archeology is the estimation of the date of an archeological event that is associated with a set of contemporaneous artifacts. This problem comes up frequently when constructing chronologies in archeology. Each artifact of the event is dated, and so we collect a Date: December 10, 2015 1 CNRS IRAMAT-CRPAA, Université Bordeaux-Montaigne and Géosciences-Rennes, Université Rennes 1.

2 Laboratoire de mathématiques Jean Leray, Unversité Nantes.

1 sample of observations provided by one or several dating methods. This leads to the issue of combining measurements in order to define the event date.

The simplest statistical model is defined as n independent measurements M 1 , . . . , M n assumed to have the same unknown mean µ. The parameter µ is in turn related to calendar time via so-called calibration curve. We denote as s 2 i the experimental variance on each measurement M i , that is supposed known and evaluated by the laboratory during the measurement process. To summarize, we write

(1)

M i = µ + s i i , ∀ i = 1, ..., n
where 1 , ..., n are n independent identically distributed Gaussian random variables with zero mean and variance 1. The maximum likelihood estimate of µ [START_REF] Ward | Procedures for comparing and combining radiocarbon age determinations: a critique[END_REF] is given by

(2) μn = n i=1 M i s 2 i n i=1 1 s 2 i .
The measurements M i with high variance s 2 i are penalized, and thus they contribute less to the estimation of the mean µ. In [START_REF] Galbraith | Optical dating of single and multiple grains of quartz from Jinmium rock shelter, Northern Australia: Part i, experimental design and statistical models[END_REF], this model is named the "common age model" and is applied for combining paleodoses in the luminescence dating method. This method is also implemented for combining 14C ages (see Bronk Ramsey, 2009a).

To take into account individual effects, we can add a random effect on the parameter µ as follows:

M i = µ i + s i i (3) µ i = µ + σλ i ,
where λ 1 , ..., λ n , 1 , ..., n are 2n independent and identically distributed Gaussian random variables with zero mean and variance 1. Error σλ i represents the uncertainty between the measurements and the event µ due to sampling (representativeness) problems of unknown origin which are not related to the measurement process in the laboratory. This model is named the "central Age model" in [START_REF] Galbraith | Optical dating of single and multiple grains of quartz from Jinmium rock shelter, Northern Australia: Part i, experimental design and statistical models[END_REF]. An explicit form of the likelihood estimate is not available. Alternatively, a Bayesian approach can be adopted to estimate the parameters µ and σ 2 of this hierarchical model. It is thus necessary to choose a prior distribution on (µ, σ 2 ). This choice is discussed for instance in [START_REF] Congdon | Applied Bayesian hierarchical methods[END_REF] or [START_REF] Spiegelhalter | BAyesian Approaches to Clinical Trials and Health-Care Evaluation[END_REF] in the particular case of meta-analysis.

In archeology, dating laboratories might provide measurements for • a 14C age, • a paleodose measurement (TL/OSL),

• an inclination, a declination or an intensity of the geomagnetic field (AM).

The laboratory measurements are then converted into calendar dates using a calibration curve (see Section 2). This step must be added to the combination model. We thus propose to extend the model with the random effect defined in (3) by adding a calibration step, and including individual effects on the variance by replacing σ by σ i .

The inclusion of individual effects is motivated by the fact that each measurement can be affected by irreducible errors [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] which can come from different sources such as:

(1) The way of ensuring that the samples studied can realistically provide results for the events that we wish to characterize (measurement or date) (2) The care taken in sampling in the field, (3) The care taken in sample handling and preparation in the laboratory, (4) Other non-controllable random factors that can appear during the process.

The rest of the paper is organized in three sections. Section 2 is devoted to the Bayesian calibration of measurements and the combination of dates. In Section 3, we describe our hierarchical Bayesian model for combining measurements in the context of dating problems and define an "event model" to estimate a date for contemporary artifacts. Section 4 provides an application to wiggle-match dating of tree-ring sequences.

Remark 1. We define each Bayesian model using a directed acyclic graph (DAG). This graph describes the dependencies in the joint distribution of the probabilistic model. Each random variable of the model (that is an observation or a parameter) appears as a node in the graph. Any node is conditionally independent of its non-descendents given its parents. Nodes represented by circles correspond to the random variables of the model. With the color of the circles, we distinguish between observations (red) and parameters (blue). The green squares indicate exogenous variables.

Calibration and linked models

In this section we describe two standard Bayesian models:

• the individual calibration of a measurement provided by a dating laboratory, • the combination of measurements performed on the same artifact using the same dating method.

2.1. Calibration. The simplest problem is the calibration of an individual measurement M in calendar time t. The measure M has a Gaussian distribution with mean µ and known variance s 2 . Following a hierarchical model, the variable µ has a Gaussian distribution with mean g(t) and variance σ 2 g (t), where g is a function called "calibration curve" linking the measurement to calendar time t. Every dating method has its own calibration curve. For instance, the curve IntCal13 [START_REF] Reimer | Intcal13 and marine13 radiocarbon age calibration curves 0-50,000 years cal bp[END_REF] converts a radiocarbon Age in calendar date for the Northern Hemisphere atmospheric samples. The calibration of 14C ages is detailed by Bronk Ramsey (2009a); [START_REF] Buck | The BAyesian Approach to Interpreting Archaeological Data[END_REF]; [START_REF] Litton | The Bayesian approach to the interpretation of archaeological data[END_REF]. In archaeomagnetism the secular variation curves of the Earth's magnetic field are used to convert an inclination, declination or intensity measurement to a calendar date (see [START_REF] Lanos | Bayesian inference of calibration curves: Application to archaeomagnetism[END_REF]). In the case of an age M provided by luminescence (TL or OSL), the calibration function is merely a linear function g(t) = (t 0 -t) where t 0 is the determination year of the age measurement M in the laboratory.

The construction of the calibration curve depends on the dating method. In the context of OSL dating the calibration curve comes from a physical equation, and so σ g is assumed to be null. In radiocarbon or archaeomagnetic dating methods, the curve comes from an estimation step (e.g. a non linear regression on a set of dated reference measurements). In this case σ g corresponds to the errors due to the estimation method. Note that these curves are regularly updated according to the improvement of the reference data.

The Bayesian model for 'calibration' or conversion of a measurement M to calendar date t is the following. The model on the measurement M is given by

M = µ + s ,
where µ is the true value of the measurement and s 2 is the known variance. The distribution of is the standard Gaussian distribution.

The calibration step converts µ to calendar date t, using the relation

µ = g(t) + σ g (t)ρ,
where both functions g and σ g are assumed known, and where ρ is a standard Gaussian random variable.

Let T be the time range of interest. This interval is fixed from prior information on the historical period to which the event belongs.

We choose the uniform distribution on T as prior distribution on the parameter t,

(4) p(t) ∝ 1 T (t).
If we integrate the posterior distribution of (t, µ) with respect to latent variable µ, we get the posterior distribution of t (up to a multiplicative constant):

(5)

p(t|M ) ∝ 1 S(t) exp -1 2S 2 (t) (M -g(t)) 2 1 T (t),
where ( 6) 2.2. Calibration from multiple measurements. We observe m independent measurements M k performed on the same object. We assume that all the measurements can be calibrated with the same curve (g, σ 2 g ). This appears, for example, in the case of radiocarbon dating when the same object is analyzed by different laboratories. In this case, the measurements have to be combined before calibration (see [START_REF] Ward | Procedures for comparing and combining radiocarbon age determinations: a critique[END_REF] and Bronk Ramsey (2009a)).

S(t) 2 = s 2 + σ 2 g (t).
The Bayesian model (Fig. 4) is defined by the following distributions:

M i |µ ∼ N (µ , s 2 i ) ∀ i = 1, ..., m µ|t ∼ N (g(t) , σ 2 g (t)) t ∼ Unif(T ),
where T is the time range of interest. This model is called R-Combine in Bronk Ramsey (2009a). One can easily get the marginal posterior density of t, which is given by

p(t|M 1 , ..., M m ) ∝ 1 Sm (t) exp -1 2 Sm (t) 2 (M m -g(t)) 2 1 T (t),
where

Sm (t) 2 = s 2 m + σ 2 g (t) with 1 s 2 m = m k=1 1 s 2 k ,
and where

M m = m k=1 M k s 2 k / m k=1 1 s 2 k .
The R-combine model is equivalent to the individual calibration of the observation (M m , s 2 m ) using the process described in Section 2.1.

t µ i = 1 to m s 2 i M i

Event model

In the general case we observe n measurements M i such that each measurement provides a dating through a calibration step defined by a calibration function g i and its error σ gi . The R-combine model is no longer valid in this case because it requires a common calibration curve (i.e. g i = g for all i = 1, ..., n). Our idea is to adapt the Bayesian combination of measurements to estimate the date of an archeological event from the individual dating of contemporaneous artifacts.

3.1. The model. We describe the so-called Event model to estimate a date θ from n measurements M i provided by different dating methods. We assume that each measurement M i is related to an individual date t i through a calibration curve g i . Here this curve is assumed known with some known uncertainty. The main assumption in our Event model is the contemporaneity of the dates t i , i = 1..., n with the event date θ. We assume that θ ∈ T where T is the bounded interval corresponding to the range of study.

In this context the model with random effect given by (3) can be rewritten as follows

M i = µ i + s i i , µ i = g i (t i ) + σ gi (t i )ρ i , t i = θ + σ i λ i , (7)
where ( 1 , ... n , ρ 1 , ..., ρ n , λ 1 , ..., λ n ) are independent and identically distributed Gaussian random variables with zero mean and variance 1. The random variables (λ i ) i and ( i ) i are interpreted as follows :

• σ i λ i represents the irreducible error between t i and θ due to sampling problems external to the laboratory ( Section 1), • s i i + σ gi (t i )ρ i represents the experimental error provided by the laboratory and the calibration step. According to the DAG defining the event model (Fig. 5), the joint distribution of the probabilistic model can be written in the form ( 8)

θ i = 1 to n σ 2 i t i µ i s 2 i M i
p(M 1 , ..., M n , µ 1 , ..., µ n , t 1 , ..., t n , σ 2 1 , ..., σ 2 n , θ) = p(θ) n i=1 p(M i |µ i )p(µ i |t i )p(t i |σ 2 i , θ)p(σ 2 i ),
where the conditional distributions that appear in the decomposition are given by

M i |µ i ∼ N (µ i , s 2 i ) µ i |t i ∼ N (g i (t i ) , σ 2 gi (t i )) t i |σ 2 i , θ ∼ N (θ , σ 2 i ) (9) σ 2 i ∼ Shrink(s 2 0 ) (10) θ ∼ Unif(T ). ( 11
)
The density of the uniform shrinkage distribution with parameter s 2 0 (denoted Shrink(s 2 0 )), that appears in (10), is given by

(12) p(σ 2 i ) = s 2 0 (s 2 0 + σ 2 i ) 2 .
Note that this form of density implies that the random variable s 2 0 /(σ 2 +s 2 0 ) is uniformly distributed on [0, 1]. This distribution is rather diffuse in the sense that the variance and mean are infinite. So it can be considered weakly informative. The properties of the uniform shrinkage prior and the choice of s 2 0 are discussed in [START_REF] Daniels | A prior for the variance in hierarchical models[END_REF]. In the particular case of the central age model defined in (3), [START_REF] Spiegelhalter | BAyesian Approaches to Clinical Trials and Health-Care Evaluation[END_REF] suggests the following choice of parameter:

1 s 2 0 = 1 n n i=1 1 s 2 i .
Parameter s 2 0 quantifies the magnitude of error on the measurements. As s 2 0 is the median of the uniform shrinkage prior, this choice ensures that the parameter σ 2 i has the same prior probability to be smaller or larger than s 2 0 . Therefore we do not favor the measurement errors with respect to the error between t i and θ due to sampling problems.

With calibrated measurements, the parameter s 2 0 cannot be calculated directly from the variances s 2 i . Indeed the measurement errors s i are not necessarily homogeneous units (see for instance the case of archaeomagnetism combined with 14C). Moreover they do not contain the information on the variance of the combined variables. To adapt this strategy we have to estimate the variance of the dates t i . We propose to estimate these variances after the step of individual calibration. We proceed in the following way, for each i = 1, ..., n:

(1) An individual calibration step is performed for each measurement M i (a) compute the posterior distribution of t i (using ( 5)) (b) approximate the posterior variance w 2 i = var(t i |M i ) (2) Take as shrinkage parameter s 0 :

1 s 2 0 = 1 n n i=1 1 w 2 i .
Remark 2. For the hyperparameters σ 2 i in the third stage of the hierarchical model, the choice of a diffuse noninformative prior may be problematic as improper priors may induce improper posteriors. This problem appears in our model, for instance by taking p(σ i ) ∝ 1/σ 2 i , that is the classical non informative prior (called the Jeffreys prior) for Gaussian observations. The choice of diffuse priors in such models raises particular issues, discussed for instance in [START_REF] Congdon | Applied Bayesian hierarchical methods[END_REF]. It is always possible to choose a proper probability distribution which approximates the diffuse prior, for instance

• a density proportional to the Jeffreys prior truncated on an closed interval [a, b] ⊂]0, ∞[, • an inverse Gamma distribution with small scale and shape parameters.

However, such priors may cause identifiability problems as the posteriors are close to being empirically improper. Moreover the inference of the parameter of interest is sensitive to the choice of parameters in the proper prior distribution. This is a major drawback for small samples.

MCMC algorithm.

The posterior distribution of the parameter of interest θ can not be obtained explicitly. It is necessary to implement a computational method to approximate the posterior distribution, its quantiles, the Bayes estimates and the highest posterior density (HPD) regions. We adopt a MCMC (Markov Chain Monte Carlo) algorithm known as the Metropoliswithin-Gibbs strategy because the full conditionals cannot be simulated by standard random generators. For each parameter, the full conditional distribution is proportional to (8).

1. For the parameter of interest θ we identify a truncated Gaussian distribution on the period of interest T . Such a distribution can be simulated using rejection sampling. Different choices of proposal distribution are possible, such as the Gaussian distribution or the truncated Laplace (double exponential) distribution on T . It is also possible to use an adaptive random walk Metropolis-Hastings (MH) with a Gaussian proposal distribution.

2. The density of the full conditional distribution of σ 2 i is explicitly computable up to an unknown multiplicative constant, but it is not a standard distribution. Therefore it is simulated using an adaptive random walk MH with a Gaussian proposal distribution on the variable log(σ 2 i ). 3. The full conditional distribution of t i , i = 1, ..., n, is proportional to

1 S i (t i ) exp -1 2S 2 i (t i ) (M i -g i (t i ) 2 exp -1 2σ 2 i (t i -θ) 2 ) ,
where S i (t i ) is defined in (6). We can choose an adaptive random walk MH with a Gaussian proposal. However, the random walk solution is not necessarily the most efficient choice because the target distribution can be multimodal. We are frequently confronted with this problem as shown by, for instance, the results yielded by archaeomagnetic dating in Example 2. In this context, an alternative is to choose an independent Hasting-Metropolis algorithm with a proposal distribution that mimics the individual calibration density defined in ( 5). This ensures that all the possible values of t i can be visited when the calibrated date distribution is multimodal.

These algorithms are implemented in the cross-platform ChronoModel application [START_REF] Lanos | Chronomodel : Chronological Modelling of Archaeological Data using BAyesian Statistics[END_REF][START_REF] Vibet | Chronomodel user's manual[END_REF], which is free and open source software. In the examples discussed below, the graphics summarizing the numerical results are performed using ChronoModel software. We represent the marginal densities of • the parameter of interest θ (on gray background) • the individual measurements t i (on white background) • the individual standard deviations σ i For each density, the bar above the density represents the shorter 95% posterior probability interval. The vertical lines, delimiting the colored area under the density curve, indicate the endpoints of the 95% HPD region. For a unimodal posterior density, the 95% posterior probability regions coincide.

Different graphical tools are implemented to assess the convergence of the MCMC in ChronoModel software (e.g. the autocorrelation functions, the acceptance rate of Metropolis-Hasting algorithms). The user can adjust the length of burn-in, the maximum number of iterations for adaptation and for acquisition, and the thinning rate. Moreover multiple chains can be simulated and exported for further analysis using, for instance, the R package coda.

In term of computational performance, among the examples given in this paper, the longest computing time is achieved for Ex. 2, Indeed 10 6 iterations are required to well estimate the multimodal posterior distribution of the Event. The simulation is performed on a 2 GHz Intel Core i7 by using a single core. The computational time is less than 1 minute including the post-processing i.e. estimation of HPD regions, credibility regions, posterior densities, etc from simulated Markov chain, and displaying results.

Examples.

Example 1. Mont-Saint-Michel (Normandy, France)

An archeology building study has been conducted at the Mont-Saint-Michel site in Normandy, France [START_REF] Blain | TL-dating applied to building archaeology: The case of the medieval church Notre-Dame-Sous-Terre (Mont-Saint-Michel, France)[END_REF][START_REF] Sapin | Archeologie du bati et archeometrie au Mont-Saint-Michel, nouvelles approches de Notre-Dame-Sous-Terre[END_REF]. The aim of this study was to estimate the dates of different building states of the Carolingian church Notre-Dame-sous-Terre (NDST) located inside the Mont-Saint-Michel abbey. Therefore it is reasonable to fix the time range T equal to [0 , 2000], which widely includes the Carolingian period. In what follows, we model building state number 1 which is associated with 14C analysis on 3 charcoals, and luminescence analyses of 8 bricks.

Figures 6 and7 summarize the estimation results obtained by the event model. The 95% HPD interval for the date of building state 1 is [892 ; 993]. This result is more precise than the dating obtained in [START_REF] Blain | TL-dating applied to building archaeology: The case of the medieval church Notre-Dame-Sous-Terre (Mont-Saint-Michel, France)[END_REF], which is based solely on luminescence analysis. This result confirms the value of combining different dating methods.

Figure 6 provides a comparison between the individual calibrations of each measurement and the posterior densities of the dates t i . Clearly the event model improves the precision of the dating of each artifact included in the event. The graphical comparison of the posterior densities of t i confirms our assumption that the bricks and charcoals are contemporary.

Posterior densities for standard deviations σ i fall within the small range [0, 200] at 95% (Fig. 7). These results confirm that the dates t i are contemporary Example 2. Cuers (Provence, France), medieval or modern lime kiln Archeological excavation was carried out to estimate the last firing date of a lime kiln in Cuers (Provence, France) [START_REF] Vaschalde | La datation des structures de cuisson: intégration de l'archéomagnétisme et du radiocarbone, apports de l'anthracologie[END_REF]. We choose the time range T equal to [0 , 2000]. Indeed no historical or archeological information is available about this site between the Roman period and the 20th century.

The studied artifacts are pieces of baked clay from the kiln wall submitted for AM analysis and two charcoals submitted for 14C analysis. The obtained measurements are three archaeomagnetic parameters (inclination, declination and intensity) and two radiocarbon ages. Figure 8 gives the results obtained by the event model. As in the previous example, posterior densities t i (in red) are much more precise than the individual calibrated densities, in the sense that the number of local modes is reduced. The Event model gives a 95% HPD region for the last firing date which is the union of two disjoint intervals [1388, 1548] and [1556, 1616]. In this case the dating is indeterminate in the sense that two periods remain possible for the last firing. Consequently, additional measurements or prior information are required to decide between these two solutions.

Figure 9 gives the posterior densities of standard deviations σ i . It clearly shows that this parameter takes higher values in the case of archaeomagnetic dating. This indicates a lack of consistency between the dates within the event. The archaeomagnetic dates are penalized in the estimation of the event. This is explained by the fact that radiocarbon dating provides more accurate dates after the calibration step. Posterior densities obtained for standard deviations σ i .

We consider the example studied in Bronk [START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF], in order to compare the robustness to outliers of the event model with the R-combine model with the outlier implemented in OxCal software.

From eleven radiocarbon dates, the aim is to determine the date of the context X in Tell Qasile, an archeological site in Tel Aviv. Since the site belongs to the protohistoric period, we fix the time range T equal to [-2000 , 0].

In OxCal, it is necessary to fix the prior probability that each measurement is an outlier. The recommended value is 5% when no pertinent information is available. The 95% HPD region for the date of the context X yielded by OxCal is [-1054; -970] ∪ [-962; -934].

The event model gives a uni-modal posterior density and 95% HPD interval equal to [-1050; -951] (Fig. 10). This is very similar to the two merged OxCal intervals.

Posterior densities for standard deviations σ i are very similar (Fig. 10) except for samples QS2 and QS6 which show higher posterior values, and thus appear to be outliers. The same two outliers are also detected by OxCal.

In conclusion, the event model is robust to outliers and has the advantage that it does not require prior information about the outliers. The ability of the individual variance σ 2 i to take large values automatically penalizes an outlier.

Wiggle-matching

The "wiggle-matching" model combines radiocarbon dating and dendrochronology. Radiocarbon dating is carried out on tree-ring samples separated by a known number of tree-rings. This gives prior information on the calibrated dates t i . The number of years between these dates is then known, [START_REF] Manning | 14C record and wiggle-match placement for the Anatolian (Gordion Area) Juniper treering chronology 1729 to 751 cal BC, and typical Aegean/Anatolian (growing season related) regional 14C offset assessment[END_REF][START_REF] Galimberti | Wiggle-match dating of tree-ring sequences[END_REF][START_REF] Litton | A Bayesian approach to wiggle-matching[END_REF].

In this case, the calibrated dates t i should be shifted according to their relationship to the event date θ. So we adapt our event model to wiggle-matching as follows.

We consider that the event date θ corresponds to the date of a chosen reference tree-ring (for instance the oldest tree-ring, see Example 4).

We denote by δ i the number of years between t i and θ, which is fixed by counting the number of rings between them and assuming that one ring is equal to one year. By convention, δ i is positive (resp. negative) when t i is older (resp. younger) than θ.

Equation ( 7) can be rewritten as follows (13)

t i = θ -δ i + σ i λ i .
Remark 3. It is also possible to assume δ i unknown in order to model a counting error in the number of tree-rings. This extension of the model is implemented in chronomodel software. An application of the Event wiggle-matching model is presented here in the case of the Gordion Juniper dendrochronology (Manning and Kromer, 2011, table 1). Thirty-five samples of this sequence, from relative tree-rings centered at 776.5-1025.5 have been dated by the Heidelberg radiocarbon laboratory. The time range T is fixed equal to [-5000 , 0], the lack of prior information leads to the choice of a very large interval.

θ i = 1 to n δ i σ 2 i t i µ i s 2 i M i
The samples are calibrated with IntCal04.14c curve [START_REF] Reimer | Intcal04 terrestrial radiocarbon age calibration, 0-26 cal kyr bp[END_REF] and are separated with known gaps taking values between 1 and 19 years over a range of 249 years. δ i is the gap between the dated tree-ring and reference. Following [START_REF] Manning | Radiocarbon dating archaeological samples in the eastern mediterranean, 1730 to 1480 bc: further exploring the atmospheric radiocarbon calibration record and the archaeological implications[END_REF], we fix as reference, the oldest tree-ring i.e. δ 1 = 0.

The model considered in [START_REF] Manning | Radiocarbon dating archaeological samples in the eastern mediterranean, 1730 to 1480 bc: further exploring the atmospheric radiocarbon calibration record and the archaeological implications[END_REF] is implemented in OxCal, by using the D-Sequence and Gap functions with the outlier model. It gives a 95% HPD interval equal to [-1734, -1724].

The Event model gives a 95% HPD interval for θ equal to [-1744, -1719], which is less precise than the estimate yielded by OxCal. The difference in the 95% HPD region is due to the event structure, which increases the variance of the posterior distribution. However our approach brings robustness. As an illustration, we artificially modify the dataset by adding outliers with the same experimental variances on the 14C ages (see Table 1). Figure 12 gives the estimation obtained by the event model for:

• the original Gordion dataset analyzed in [START_REF] Manning | Radiocarbon dating archaeological samples in the eastern mediterranean, 1730 to 1480 bc: further exploring the atmospheric radiocarbon calibration record and the archaeological implications[END_REF],

• the same dataset contaminated with 11 outliers. Manning and Kromer (2011) notations 20144, 27605 3100 20137, 25792, 20155 3600 20157, 20141, 20147, 25793, 20153, 27612 3700 Figure 12 shows that the posterior distribution of the event date is very weakly sensitive to the presence of these outliers. Indeed, on the contaminated sample, the event model gives an 95% HPD interval for θ equal to [-1749 ; -1717]. The two results are quite similar; the event date is not affected by the presence of the outliers. The event model considered in [START_REF] Manning | Radiocarbon dating archaeological samples in the eastern mediterranean, 1730 to 1480 bc: further exploring the atmospheric radiocarbon calibration record and the archaeological implications[END_REF] cannot be estimated using OxCal software because of poor agreement indices. This alert message indicates in particular that the convergence of the MCMC is not guaranteed.

Figure 12. Gordion juniper dendrochronology (Central Anatolia). Posterior densities of the dates t i (solid line and white background), the shifted dates t i + δ i (dotted line and white background) and the Event date θ (gray background) obtained on original Gordion dataset (blue) and on a contaminated version with 11 outliers defined in Table 1 (in red).

Perspective

A combination procedure to estimate the date of an archeological event from a set of contemporaneous artifacts using a hierarchical Bayesian approach is proposed. The event model defines a fundamental chronological entity that can be related to other entities in various ways to produce complex chronological models. We can then construct chronologies of events which take into account new archeological information such as ordering constraints based on stratigraphy, phasing with succession or with a hiatus between phases, or with a duration constraint on the phase, etc. These features are implemented in the ChronoModel software.

Figure 1 .

 1 Figure 1. Conversion of a 14C age (A = 1000 ± 30 BP) to a calendar date via the calibration curve IntCal13.Figures 1, 2 and 3 illustrate the calibration process for different calibration curves. In many cases, the individual calibration provides an estimation of the date with a high uncertainty. The posterior distribution is often multi-modal or unimodal with large variance. This motivates the archeologist to repeat measurements on different contemporaneous artifacts using one or several dating methods.

Figure 2 .

 2 Figure 2. Conversion of a TL age (A = 900 ± 50 ) to a calendar date through the linear transformation g(t) = (2000 -t).

Figure 3 .

 3 Figure 3. Conversion of an inclination measurement (I = 59 ± 1) to a calendar date via the calibration curve of archaeomagnetic field in France (Paris) over the last two millennia.

Figure 4 .

 4 Figure 4. DAG of the R-Combine model.

Figure 5 .

 5 Figure 5. DAG for the hierarchical Event model applied to dating combination with calibration.

Figure 6 .

 6 Figure 6. Mont-Saint-Michel, Ex. 1. Posterior density of the event date θ (gray background). Posterior densities of t i (blue line & white background) and individual posterior calibrated densities (gray line & white background) obtained for TL dates and for 14 C dates.

Figure 7 .

 7 Figure 7. Mont-Saint-Michel, Ex. 1 (Cont.). Posterior densities obtained for standard deviations σ i .

Figure 8 .

 8 Figure 8. Cuers, Ex. 2. Posterior density of the event date θ (gray background). Posterior densities of t i (red line & white background) and individual posterior calibrated densities (gray line & white background) obtained for archaeomagnetic and 14 C dates.

Figure 9 .

 9 Figure 9. Cuers, Ex 2 (Cont.) Posterior densities obtained for standard deviations σ i .

Example 3 .

 3 Figure 10. Tell Qasile, Ex. 3. [left] Posterior density for Event date θ (gray background). Posterior densities t i for 14 C dates (blue line and white background) The individual posterior calibrated densities are superimposed in gray. [right] Posterior densities obtained for standard deviations σ i .

Figure 11 .

 11 Figure 11. DAG for the wiggle matching model.

Example 4 .

 4 Gordion Juniper dendrochronology (Central Anatolia)

Table 1 .

 1 Gordion Juniper dendrochronology (Central Anatolia).

	Contaminated

Acknowledgement

The authors are extremely grateful to Thomas S. Dye for his careful reading and his insightful comments that helped to improve the exposition in the manuscript. This project is supported by the grant ANR-11-MONU-007 ChronoModel.


Philippe Lanos
email: philippe.lanos@univ-rennes1.fr

Anne Philippe

EVENT MODEL: A ROBUST BAYESIAN TOOL FOR CHRONOLOGICAL MODELING

Keywords: Event model, hierarchical Bayesian statistics, individual errors, outlier penalization, robustness, MCMC methods, Chronomodel software

We propose a new modeling approach for combining dates through the Event model by using hierarchical Bayesian statistics. The Event model aims to estimate the date of a context (unit of stratification) from individual dates assumed to be contemporaneous and which are affected by errors of different types: laboratory and calibration curve errors and also irreducible errors related to contaminations, taphonomic disturbances, etc, hence the possible presence of outliers. The Event model has a hierarchical structure which makes it possible to distinguish between date of an Event and dates of the artifacts involved. Prior information on the individual irreducible errors is introduced using a uniform shrinkage density with minimal assumptions about Bayesian parameters. The model is extended in the case of stratigraphic sequences which involve several Events with temporal order constraints (relative dating). Calculations are based on MCMC numerical techniques and can be performed using the ChronoModel software which is freeware, open source and cross-platform. This modeling provides a very simple way to automatically penalize outlying data without having to remove them from the dataset. This approach is compared to alternative approaches implemented in Oxcal or BCal software: we show that the Event model is more robust but generally yields less precise credibility intervals. Mathematical formulations are explained in detail and comparisons are done thanks to synthetic examples. Three application examples are shown: the radiocarbon dating of the shroud of Turin, the dating of a medieval potter's kiln in Lezoux (Auvergne, France) by using radiocarbon, archaeomagnetism and thermoluminescence, and the OSL dating of the Shi'bat Dihya 1 sequence in Wadi Surdud middle paleolithic complex (western Yemen).

Introduction

Bayesian chronological modeling appears as an important issue in archeology and palaeoenvironmental sciences. It is now standard practice especially in England [START_REF] Bayliss | Rolling out revolution: Using radiocarbon dating in archaeology[END_REF] where the methodology has been developed since the 1990s. Internationally, Bayesian chronological modeling has been adopted more slowly, but it is now beginning to be the method of choice for the interpretation of radiocarbon dates, at least in the specialist literature. This can be illustrated, for example, by statistical methods used for the interpretation of radiocarbon dates in papers published in the journal "Radiocarbon" arising from the regular series of "Radiocarbon" and "14C in Archaeology" conferences [START_REF] Bayliss | Quality in bayesian chronological models in archaeology[END_REF]. Most applications are undertaken using the flexible software packages, BCal [START_REF] Buck | BCal : an on-line Bayesian radiocarbon calibration tool[END_REF], Datelab [START_REF] Nicholls | New radiocarbon calibration software[END_REF] and OxCal [START_REF] Bronk Ramsey | Radiocarbon calibration and analysis of stratigraphy : the OxCal program[END_REF][START_REF] Bronk Ramsey | Probability and dating[END_REF][START_REF] Bronk Ramsey | Development of the radiocarbon calibration program OxCal[END_REF], 2008, 2009a,b;Bronk Ramsey et al., 2001, 2010;[START_REF] Bronk Ramsey | Recent and planned developments of the program OxCal[END_REF]. Calibrated radiocarbon dates, or other date estimates on the calendar scale, are combined with prior archeological information of various kinds to produce a combined chronology that should be more reliable than its individual components. Chronological models are thus interpretative constructions. The aim is to estimate the date of events or the chronological timing of phases. In Bayesian chronological modeling, a phase is a group of contexts (units of stratification) of similar age. This is an interpretative construct formulated with both stratigraphic and non-stratigraphic information. Thus we can distinguish a stratigraphic phase and a chronological phase (period) [START_REF] Dye | Archaeological sequence diagrams and bayesian chronological models[END_REF]. The phase is characterized by a beginning, an end and therefore a duration. We need to date depositional or interfacial contexts belonging to the phase in order to be able to estimate these temporal parameters.

Contrary to the phase which suggests a duration, we define the Event as the date of an archeological context determined from a collection of contemporaneous artifacts, which are dated by using different dating methods such as radiocarbon, thermoluminescence, archaeomagnetism, etc [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. The dating process provides measurements which are converted into calendar dates using calibration reference curves. All the dates in an Event are assumed contemporary but they can be disturbed by errors of unknown origin which can arise from different type of sources [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF]:

(1) the way of ensuring that the samples studied can realistically provide results for the events that we wish to characterize (measurement or date), (2) the care in sampling in the field, the care in sample handling, preparation and measure in the laboratory, (3) other non-controllable random factors that can appear during the process.

We propose a model of chronology called Event in the sense that:

• dates of assumed contemporaneous artifacts are embedded in one Event,

• in stratigraphic sequences, constraints of chronological order are defined between Events, not between individual dates of the artifacts.

This paper concentrates on the Event model, considered alone or incorporated in the stratigraphic sequence. Its statistical properties are described in Section 2. Alternative models developed by [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF], [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF][START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF][START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF] are presented in Section 3 and compared with our chronological model. Application examples are discussed in section 4.

The Event model

The Event model we propose for combining dates aims to be robust to outliers in the sense that the estimates are not sensitive to the presence of outliers. Consequently, our approach does not model the outliers, i.e. we do not estimate the posterior probability that a date is an outlier (See Section 3 for such an approach). Generally such a robust approach leads to a loss of accuracy in posterior estimates (for instance we obtain longer HPD intervals). Outlier modeling can provide more accurate results, but it often requires two (maybe more) estimations of the model. Indeed the outliers are identified after a first estimation and thus discarded from the dataset. Then the final model is estimated again from the new dataset. The Event model which is based on the choice of robustness, avoids this two-step procedure. It is described in the following section.

2.1. The statistical model. The Event θ is estimated from n independent measurements (observations) M i . Each measurement M i , obtained from a specific dating technique, is related to an individual date t i through a calibration curve g i and its error σ gi (see Section 2 in [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. Here this curve is supposedly known with some known uncertainty. If we assume that all the measurements can be calibrated with a common calibration curve (i.e. g i = g for all i = 1, ..., n), for example when the same object is analyzed by different laboratories, the measurements can be combined beforehand according to the R-combine model [START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF]. On the contrary, if several calibration curves are involved, the R-combine model is no longer valid and can be replaced by the Event model. The main assumption we make in the Event model is the contemporaneity of the dates t i , i = 1..., n with the event date θ. Because of error sources of unknown origin, there exists an 'overdispersion' of the dates with respect to θ. This overdispersion which corresponds to irreducible errors [START_REF] Niu | The Bayesian approach to radiocarbon calibration curve estimation: the intcal13, marine 13, and shcal13 methodologies[END_REF] is modeled by individual errors σ i .

In this context the model with random effect can be written as follows

M i = µ i + s i i , µ i = g i (t i ) + σ gi (t i )ρ i , t i = θ + σ i λ i , (1)
where ( 1 , ... n , ρ 1 , ..., ρ n , λ 1 , ..., λ n ) are independent and identically Gaussian distributed random variables with zero mean and variance 1. The random variables (λ i ) i and ( i ) i are independent and satisfy the following properties:

• σ i λ i represents the error between t i and θ due to sampling problems external to the laboratory (see error sources in Section 1), • s i i + σ gi (t i )ρ i represents the experimental error provided by the laboratory and the calibration error. The joint distribution of the probabilistic model can be written according to a Bayesian hierarchical structure:

p(M 1 , ..., M n , µ 1 , ..., µ n , t 1 , ..., t n , σ 2 1 , ..., σ 2 n , θ) = p(θ) n i=1 p(M i |µ i )p(µ i |t i )p(t i |σ 2 i , θ)p(σ 2 i ), (2)
where the conditional distributions that appear in the decomposition are given by:

M i |µ i ∼ N (µ i , s 2 i ), µ i |t i ∼ N (g i (t i ), σ 2 gi (t i )), t i |σ 2 i , θ ∼ N (θ, σ 2 i ), θ ∼ Uniform(T ), σ 2 i ∼ Shrink(s 2 0 ). (3)
The parameter of interest θ is assumed to be uniformly distributed on an interval T = [θ m , θ M ]. This interval, called "study period", is fixed by the user based on historical or archeological evidences. This appears as an important a priori temporal information. Note that we do not assume that the same information is available on the dates t i . Consequently their support is the set of real numbers R.

The uniform shrinkage distribution for σ 2 i , denoted Shrink(s 2 0 ), admits as density

p(σ 2 i ) = s 2 0 (s 2 0 + σ 2 i ) 2 , (4)
where the parameter s 2 0 must be fixed. The idea to be assumed is that, a priori, σ 2 i should be of the same order of magnitude as the variances of the latent variables t i [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. One evaluates this magnitude through the n individual measurement calibrations. More precisely, the parameter s 2 0 is chosen as follows for each i = 1, ..., n:

1 ) An individual calibration step is carried out for each measurement M i . From the posterior distribution of the individual calibration on the interval for which the calibration curve is defined, we approximate its posterior variance var(t i |M i ) by w i the standard Monte Carlo approximation. 2 ) Hence the shrinkage parameter s 0 is given by:

1 s 2 0 = 1 n n i=1 1 w 2 i .
(5) Remark 1. About numerical issues: the posterior distributions of the parameter of interest θ, the dates t i and the individual variances σ 2 i cannot be obtained explicitly. It is necessary to implement a computational method to approximate the posterior distributions, its quantiles, the Bayes estimates and the HPD regions. We adopt an MCMC algorithm and more precisely the Metropolis-within-Gibbs strategy. The algorithm is implemented in the cross-plateform ChronoModel application, which is free and open source software.

From an archeological point of view, it is important to note that the Event θ corresponds to the date of a context determined from the dating of assumed contemporaneous artifacts. The posterior date t i of an artifact is the calibrated date if considered alone. When integrated in an Event model within a stratigraphic sequence, this posterior date t i corresponds to an updating of the calibrated date by accounting for the once-whole model constraints. Thus, this conceptual distinction between Event (date of a context) and the dating (date of an artifact or an interface) leads to a new approach of the chronology. For instance, let us consider several radiocarbon datings from charcoals retrieved in a pyroclastic flow emitted by a volcano. The deposit of this flow is short in time: its Event (date of deposit) is then estimated from the charcoal dates. Some of these dates can be disturbed for many reasons: old wood effects, subsequent contamination due to erosion, etc., but it remains difficult to reject them because we are unable to decide (without looking at the datings !) which of them are really unreliable. In other words, the date t i of the artifact is good (it is an Event in itself), the Event of deposit θ is good, but some of the dates are not well connected to the Event without knowing why. The Event model is constructed to manage such a situation in a robust way so that the event remains insensitive to these "unknown" outliers. Although the probability of a date to be an outlier is not estimated, it is however possible to assess the outlying character of this date looking at the posterior density of the standard deviation σ i .

2.2. Event model pertaining to relative dating. Relative dating based on stratigraphy as defined in [START_REF] Harris | Principles of archaeological stratigraphy[END_REF] and [START_REF] Desachy | Du temps ordonné au temps quantifié : application d'outils mathématiques au modele d'analyse stratigraphique d'edward harris[END_REF][START_REF] Desachy | De la formalisation du traitement des données stratigraphiques en archéologie de terrain[END_REF], implies antero-posteriority relationships between Events. In the Bayesian framework, this temporal order information can be taken into account in the construction of the prior distribution and makes it possible to significantly improve chronometric datings (see initial articles of [START_REF] Buck | Combining archaeological and radiocarbon information : a Bayesian approach to calibration[END_REF][START_REF] Buck | Calibration of radiocarbon results pertaining to related archaeological events[END_REF][START_REF] Buck | A case study in combining radiocarbon and archaeological information : the early bronze age of st-veit-klinglberg, land salzburg, austria[END_REF][START_REF] Buck | The BAyesian Approach to Interpreting Archaeological Data[END_REF] and [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF]). In our modeling, the stratigraphic relationship is placed between the Events, no longer between the dates as in Bcal [START_REF] Buck | BCal : an on-line Bayesian radiocarbon calibration tool[END_REF] or Oxcal [START_REF] Bronk Ramsey | Radiocarbon calibration and analysis of stratigraphy : the OxCal program[END_REF][START_REF] Bronk Ramsey | Recent and planned developments of the program OxCal[END_REF].

To estimate an Event, we need to incorporate several datings in it. However, it is possible to nest only one dating per Event provided that the set of Events is constrained by temporal order.

To condense the writing, one defines an Event by the synthetic notation E j which groups measurements (observations) and variables (M j , t j , σ 2 j , θ j ), with M j = (M ji q) T for i = 1, ..., n j . Here the dot symbol means that we can take several replicated measurements M jik for one artifact. Consequently, measurements are combined beforehand using the R-Combine model described in Bronk Ramsey (2009b). For a set of r Events E j submitted to temporal order, the probabilistic model for Events chronology (see DAG in Figure 1) becomes:

p(E 1 , ..., E r ) = p(θ 1 , ..., θ r ) r j=1 nj i=1 p(M ji q|µji)p(µji|tji)p(tji|σ 2 ji , θ j )p(σ 2 ji ), (6)
where

M ji q|µji ∼ N (µ ji , s 2 ji q ) µ ji |t ji ∼ N (g ji (t ji ), σ 2 gji (t ji )) t ji |σ 2 ji , θ j ∼ N (θ j , σ 2 ji ) σ 2 ji ∼ Shrink(s 2 0j ) (θ 1 , ..., θ r ) ∼ Uniform(C) with C = S ∩ T r and • T r = [θ m , θ M ] r
the support which defines the study period,

• S = the set of r-uplets Events θ j which respect total or partial order relationships. The temporal order can be determined from a stratigraphic (or sequence) diagram. This diagram has a directed acyclic graph (DAG) structure (see Remark 2). The order relation (<) is irreflexive, asymmetrical and transitive. Consequently, stratigraphic diagrams present neither symmetry nor circularity. Note that it is not necessary to indicate transitivity relationships between Events because this indication is redundant. For each Event, it is sufficient to indicate what is beneath and what is above [START_REF] Dye | Archaeological sequence diagrams and bayesian chronological models[END_REF].

Remark 2. Each Bayesian model can be described using a directed acyclic graph (DAG). Such a graph describes the dependencies in the joint distribution of the probabilistic model. Each random variable of the model (that is an observation or a parameter) appears as a node in the graph. Any node is conditionally independent of its non-descendants given its parents. Hereafter, the circles correspond to all the random variables of the model. With the color of the circles, we distinguish between observations (red), parameters (blue) and exogenous variables (green).

It is not rare to encounter dating results which contradict the stratigraphic order: one speaks of "stratigraphic inversion". This situation often occurs when some artifact movements are provoked for example by bioturbations or establishment of backfill soils. The Event model makes it possible to manage such situations thanks to the individual variances σ 2 i which automatically penalize the dates that are inconsistent with the stratigraphic order. Applications to stratigraphic sequences are shown in examples 2 and 5.

j = 1 to r and θ j ∈ S i = 1 to n j θ j t ji µ ji M ji q σ 2 ji s 2 ji q Figure 1
. DAG for the hierarchical Event model defined in (6).

Alternative chronological models

It is important to compare the Event model properties with existing alternative models. First of all, in the following, Event parameter θ directly coincides with date parameter t i . Each dating is then representative of a singular event and thus, any irreducible error directly impacts either the measurement M i or the time t i .

Several outlier models have been implemented from the 90's: A first approach was proposed by [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] and developed thereafter by Bronk Ramsey (2009b). These models are typically based on outlier modeling. We may cite Bronk Ramsey (2009b):

outlying samples are progressively down-weighted as they are more likely to be outliers and so the results from the analysis are essentially an average between a model in which the measurement is accepted and one in which it is rejected. If we do not wish to have model averaging, but do wish to use outlier analysis solely for outlier detection, we should first run a model with outlier analysis, see which measurements are likely to be spurious and then run it again, without outlier analysis but with some of the spurious results removed entirely.

An alternative approach has been proposed by [START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF] for the analysis of 14C data where the associated variance is taken as the product of an unknown constant with the sum of the variance reported by the laboratory and the variance of the calibration curve (that is, an unknown error multiplier). These different models are discussed in the next three sub-sections.

3.1. Outlier model with respect to the measurement parameter. [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] approach to detect outliers basically consists in finding the posterior probability of a measurement (i.e. radiocarbon determination) being an outlier (based on a shifted model for each determination). We describe this chronological modeling following [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] , [START_REF] Buck | Bayesian tools for tephrochronology[END_REF] and the s-and r-type outlier models proposed by Bronk Ramsey (2009b). The corresponding model with random effect can be written according to a global formulation as follows (see DAG in Figure 2): M j qk = µ j q + s j qk j qk + δ j qkφj qksj qk10 u , µ j q = g(t j q) + σ g (t j q)ρj q (7) where • ( j q1, ... j qn, ρ 1 q, ..., ρ r q) are independent and identically Gaussian distributed random variables with zero mean and variance 1. • s j qk j qk represents the experimental error provided by the laboratory and σ g (t j q)ρj q the calibration error. • the prior on φ j qk is the Bernoulli distribution with parameter p j . A priori, φ j qk takes the value 1 if the measurement requires a shift and 0 otherwise. In practice p j must be chosen and the recommended values are 0.1 in [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF]; [START_REF] Buck | Bayesian tools for tephrochronology[END_REF] or 0.05 in Bronk Ramsey (2009b). • δ j qk corresponds to the shift on the measurement M j qk if it is detected as an outlier.

The prior for δ j qk is a Gaussian distribution N (0, σ 2 δ ), or a Student distribution T with ν degrees of freedom.

• u is a scale parameter to offset δ j qk. Parameter u can be fixed (for instance 0) or prior distributed as Uniform(0,4)

Remark 3. The date t j q and the Event θ j characterize the same date. The dot symbol in t j q is added to avoid any confusion with the dates t ji . This notation makes it possible to keep the same formalism in the sense that the same index position represents the same hierarchical level in the different models throughout the paper. j = 1 to r and t j ∈ S (i = q), k = 1 to m j q t j q µ j q

M j qk δ j qk φ j qk s 2 j qk
Figure 2. DAG for outlier model with respect to the measurement parameter defined in ( 7) and ( 8).

For more simplicity, the parameter u is assumed known. In this case the joint distribution of the probabilistic model can be written in the form: p(M, µ, t, δ, φ) = p(t 1 q, ..., t r q) r j=1 p(µ j q|tj q) mj q k=1 p(M j qk|µj q, δ j qk, φ j qk)p(δj qk)p(φj qk),

where the conditional distributions that appear in the decomposition are given by: M j qk|µj q ∼ N (µ j q + δ j qkφj qksj qk, s 2 j qk ) µ j q|tj q ∼ N (g j q(tj q), σ 2 gj q (t j q)) (t 1 q, ..., t r q) ∼ Uniform(C)

δ j qk ∼ N (0, σ 2 δ ) or ∼ T (ν) φ j qk ∼ Bernoulli(p j ) with C = S ∩ T r . T r = [t m , t M ]
r is the study interval and S is the set of r-uplets of time t j q which respect the total or partial order relationships between dates. This model becomes the r-type model if we replace random error δ j qkφj qksj qk by δ j qkφj qk, that is without scaling by laboratory error [START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF], page 1038).

In Section 4, we compare this approach with the Event model by using the famous example of the shroud of Turin (see Example 3) 3.2. Outlier model with respect to the laboratory variance parameter. [START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF] propose an alternative model that provides a robust approach for the analysis of 14C data, without the need to detect outliers.

An uncertainty about the laboratory-reported variance s 2 j qk is introduced by considering the product α j qs 2 j qk , where α j q is a strictly positive unknown constant. The meaning of α j q is that of an unknown "variance multiplier" applied to the laboratory-reported variance. This formulation is equivalent to a Normal model with unknown variance, but also considering the reported error s j qk as part of the data. Then, the modified statistical model for the measurement M j qk, after integration with respect to µ j q, (see DAG in Figure 3) becomes: M j qk|µj q, α j q ∼ N (µ j q, α j qs 2 j qk ) µ j q|tj q, α j q ∼ N (g j q(tj q), α j qσ 2 gj q (t j q)) (9) j = 1 to r and t j ∈ S (i = q), k = 1 to m j q t j q µ j q M j qk α j q s 2 j qk Figure 3. DAG for outlier model with respect to the experimental variance parameter defined in (9).

Using this model, the effect of outlier observations is reduced in the posterior distributions, without the need to include additional parameters or removing determinations from the data set. The resulting posterior distribution for t i has a smoother shape in comparison to the previous Normal model. The features of this modeling are somewhat similar to the Event model in the sense that it is robust to outlier observations and other causes of overdispersed data, with far fewer parameters than with [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] and Bronk Ramsey (2009b) models. However, this model does not consider individual effects because α j q is not dependent of k, contrary to the Event model where individual effects are carried out by the variances σ 2 ji . Moreover, this model depends on the choice of the inverse gamma prior with exogenous parameters a = 3 and b = 4 for the variance multiplier. [START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF] applied this outlier modeling to the shroud of Turin. In Section 4, we compare this approach with the Event model (see Example 3) 3.3. Outlier model with respect to the time parameter. This chronological modeling, called t-type outlier model, is proposed by (Bronk Ramsey, 2009b, see page 1035-1037). The model with random effect becomes (see the DAG in Figure 4): M j qk = µ j q + s j qk j qk, µ j q = g(t j q + δ j qφj q10 u ) + σ g (t j q + δ j qφj q10 u )ρ j q (10)

where parameters ( j q1, ... j qn, ρ 1 q, ..., ρ r q) , δ j q , φ j q and u are defined in the same way as in Section 3.1. When parameter u is fixed, the joint distribution of the probabilistic model can be j = 1 to r and t j ∈ S (i = q), k = 1 to m j q t j q µ j q M j qk δ j q φ j q s 2 j qk Figure 4. DAG for outlier model with respect to the time parameter defined in (10) and ( 11)

written in the form:

p(M, µ, t, δ, φ) = p(t 1 q, ..., t r q) r j=1 p(µ j q|tj q, δ j q, φ j q)p(δj q)p(φj q)

mj k=1 p(M j qk|µj q), (11)
where the conditional distributions that appear in the decomposition are given by: M j qk|µj q ∼ N (µ j q, s 2 j qk ) µ j q|tj q ∼ N (g j q(tj q + δ j qφj q), σ 2 gj q (t j q + δ j qφj q)) (t 1 q, ..., t r q) ∼ Uniform(C)

δ j q ∼ N (0, σ 2 δ ) or ∼ T (ν) φ j q ∼ Bernoulli(p j )
with C = S ∩ T r defined as in Section 3.1

This modeling can be compared to the Event model if we consider several dates i nested in the Oxcal function "Combine" with t-type outlier model. The date t j q then becomes a common date t j for measurements µ ji Consequently, ( 11) is transformed into:

p(M, µ, t, δ, φ) = p(t 1 , ..., t r ) r j=1 nj i=1 p(µ ji |t j , δ ji , φ ji )p(δ ji )p(φ ji ) mji k=1 p(M jik |µ ji ), (12)
Applying the R-Combine procedure to measurements M jik , we obtain:

p(M, µ, t, δ, φ) ∝ p(t 1 , ..., t r ) r j=1 nj i=1 p(M ji q|µji)p(µji|tj, δ ji , φ ji )p(δ ji )p(φ ji ), (13)
Particular case for comparison: If we consider linear calibration curves with constant errors, that is, by setting:

• g ji (t ji ) = t ji • σ gji (t ji ) = σ gji , and knowing that • δ j qk ∼ N (0, σ 2 δ ) • φ j qk ∼ Bernoulli(p j ), it is possible to analytically integrate (13) with respect to δ ji , φ ji , and µ ji . The posterior probability density of t is then given by:

p(t|M ) ∝ p(t 1 , ..., t r ) r j=1 nj i=1 p(M ji q|tj), (14)
where the conditional distribution of M ji . given t j is a finite mixture distribution defined by

p(M ji. |t j ) = p j 1 2π(S 2 ji q + σ 2 δ ) e 1 -2(S 2 ji q +σ 2 δ ) (Mji q-tj) 2 + (1 -p j ) 1 S ji q √ 2π e 1 -2S 2 ji q (Mji q-tj) 2 (15)
with S 2 ji q = s 2 ji q + σ 2 gji . On the other hand, posterior density of Event θ in ( 6) can be compared to density of time t j in ( 14) after integration with respect to µ ji and σ 2 ji . The posterior probability density of θ is given by:

p(θ|M ) ∝ p(θ 1 , ..., θ r ) r j=1 nj i=1 p(M ji q|θj), (16)
where

p(M ji. |θ j ) = ∞ 0 1 2π(S 2 ji q + σ 2 ji ) e 1 -2(S 2 ji q +σ 2 ji ) (Mji q-θj) 2 s 2 0j (s 2 0j + σ 2 ji ) 2 dσ 2 ji (17)
A graphical representation of the densities defined in ( 15) and ( 17) are shown in figure 5 with the following parameter values: s ji = 30 , σ gji = 10, σ δ = 10 2 and θ j = t j = 1000. The density (15) is plotted for three different values p j = 0.01, 0.05, 0.10 . We can observe that shrinkage modeling in (17) leads to a more diffuse density making it possible to better take into account the possible presence of outliers. This behaviour is illustrated in the five following examples. 15) and ( 17) with the following parameter values: s ji = 30 , σ gji = 10, σ δ = 10 2 and θ j = t j = 1000

3.4. Comparisons with synthetic data.

Example 1. Synthetic data with outlier In order to illustrate these properties, we consider a synthetic (fictitious) Event example composed of four Gaussian distributed measurements around 800, 850, 900 and 1500 with standard deviation 30, and calibrated using the simple linear transformation M = g(t) = t. The time range (study period) T is set equal to [0 , 2000], which widely includes the calibrated data. This example (Fig. 6 ) shows the effect of an outlying date (here 1500) on the event θ. The 95% HPD interval obtained with the four data is [795, 914] AD while the 95% HPD interval obtained with the three well grouped data is [791,909]. The outlying date has no influence on Event results in the sense that both intervals stay centered around the same value. However we lose some precision in the presence of the outlier.

In Figure 7, three of the posterior distributions of the standard deviations σ i (i=1,..,3) remain near zero while the standard deviation of the outlying data is high because of its outlying position. This behavior demonstrates that the Event model appears to be a robust statistics for calculating the posterior mean of the date θ with a very weak assumption on prior densities.

In Figure 6, right, the t-type outlier model is applied to the same data using the "Combine" function and the outlier model ("General",T(5),U(0,5),"t") with prior outlier probability p j = 0.10. The calculation was performed with OxCal software V4.2. The combination yields to a bimodal posterior density for t with a 95% HPD region equal to [800, 900] and [1450, 1540] AD. The posterior outlier probabilities are respectively equal to 24, 21, 24 and 80%. According to OxCal recommendations we remove the date "1500" and we rerun the estimation of the model. Then the new 95% HPD interval is [810, 890] AD, a result included in the Event interval [795,914]. This illustrates the fact that we lose in precision in order to have a robust estimation. Example 2. Synthetic data in a stratigraphic sequence with date inversion.

We consider a stratigraphic sequence composed of five Events such that:

F 5 ≤ F 4 ≤ F 3 ≤ F 2 ≤ F 1
Each Event comprises two Gaussian distributed dates, having the same standard error 30 and calibrated using the simple linear transformation M = g(t) = t:

Event(F5) Date("D51",650,30) Date("D52",600,30)

Event(F4) Date("D41",350,30) Date("D42",400,30)

Event(F3) Date("D31",700,30) Date("D32",750,30)

Event(F2) Date("D21",250,30) Date("D22",200,30)

Event(F1) C-Date("D11",800,30) C-Date("D12",850,30) The time range (study period) T is set equal to [-1000 , 2000], which widely includes the calibrated data. This example (Fig. 9, left) shows the effect of outlying dates which induce some stratigraphic inversions between Events. A posteriori results show that Event F 2 is well corrected (dates F 21 and F 22 appear as clear outliers). Conversely, Events F 4 and F 5 show ambiguous results, the posterior densities being bimodal. This result is not surprising because the dates in both Events F 4 and F 5 are earlier than the dates in Event F 3. Therefore the stratigraphic relationship does not provide information making it possible to decide which of the Events F 4 and F 5 is in stratigraphic inversion.

In Figure 9 [right], the t-type outlier model is applied to the same dataset using the "Combine" function and the outlier model ("General",T(5),U(0,5),"t") with prior outlier probability p j = 0.05. The calculation was performed with Oxcal software V4.2. The combination with the stratigraphic order yields to monomodal posterior densities for t. Events F 2 and F 4 are corrected in the same way and there is no ambiguity: only one stratigraphic scenario is retained.

However, if the prior outlier probability is modified, p j = 0.10, then some bimodal posterior densities appear for Events F 2 and F 4 , leading to a similar result to example 1. Thus, the t-type outlier model is sensitive to the choice of p j value, while the Event model does not depend on any parameter choice.

Application examples

Archeological applications of the Event model are shown in this section: in a first example, the model is applied to radiocarbon dating results from the shroud of Turin (Italy). In the second example (Lezoux site, France), the Event is determined from a set of results involving radiocarbon, thermoluminescence and archaeomagnetic dating techniques. The third example deals with the question of Events constrained within a stratigraphic sequence (Shi'bat Dihya 1 site, Yemen).

Example 3. Shroud of Turin (Italy)

Twelve radiocarbon datings have been performed on a strip cut from the shroud of Turin (Damon, 1989) and divided into three samples sent to Arizona, Oxford and Zurich AMS laboratories. Four, three and five determinations were made respectively (see page 498-499 Christen, 1994). The shroud belonging to the historical period, the prior time range T is set equal to [0 , 2000].

Figure 10 gives the posterior results with the Event model. The 95% HPD interval for the event θ is equal to [1265 ; 1315] AD. Figure 11 gives the posterior densities of standard deviations σ i . All the densities are very similar, and thus they do not clearly indicate the presence of outliers in this data set.

Comparison with outlier model with respect to the measurement parameter defined in(8) . The HPD interval with the Event model is very close to the interval [1270 ; 1310] obtained by Christen [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF] after removing the determinations A1.1 and O1.1 identified as outliers. This clearly shows that the event model can replace the R-Combine outlier model average. Comparison with outlier model with respect to the laboratory variance parameter defined in (9). According to [START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF], the 95% HPD probability interval obtained in which the organic matter in the Shroud of Turin died is [1281 ; 1302] AD. This interval is included in the one obtained by the Event model. This is explained by the fact that the Christen and Perez model considers a common parameter α for the dates of each laboratory (accordingly, the DAG in Figure 3 has to be slightly modified), while the Event model considers individual variances for all the dates, thus leading to higher robustness but less precision. The example of the medieval potter's kiln of Maison-de-Retraite-Publique site [START_REF] Mennessier-Jouannet | Convergence de la typologie de céramiques et de trois méthodes chronométriques pour la datation d'un four de potier à lezoux (puy-de-dôme)[END_REF], in Lezoux (Auvergne, France), shows how to combine data from different dating techniques using different calibration curves. Here, the Event corresponds to the last firing of the kiln and is determined from baked clay (AM and TL dating) and from charcoals (14C dating) of trees assumed to be felled at the same time as the last firing. More precisely, dating is based on 3 TL datings (CLER 202a, 202b, 203), 2 AM datings (inclination and declination) and 1 radiocarbon dating (Ly-5212) (Fig. 12). As the kiln belongs to the historical period, the prior time range T is set equal to [0 , 2000]. Posterior densities t i (in color) obtained are greatly shrunken compared to individual calibrated densities (in black line), especially for archaeomagnetic and TL dates. Event model gives a 95% HPD interval equal to [575,888] AD. Posterior densities for standard deviations σ i (Fig. 13) are much more spread than in the previous example of Turin. This comes from the multimodality of AM calibrated dates obtained with inclination (Inc) and declination (Dec). More generally, the parameter σ i takes larger values when the associated date t i is a possible outlier (see also examples in [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. In the previous examples 3 and 4, we can see that the Event model gives a very simple procedure to combine different dating results assumed to be conditionally independent to the Event date. There is no need to model the outliers nor to adjust some specific parameters for outlier detection. When the dates are coherent, posterior densities for standard deviations remain similar and concentrated near zero. In the next example, we consider a stratigraphic sequence of sandy silt deposits dated by OSL.

Example 5. Wadi Surdud paleolithic site (Yemen): OSL dating sequence. The Shi'bat Dihya 1 (SD1) site in Wadi Surdud middle paleolithic complex (western Yemen) (see [START_REF] Delagnes | Inland human settlement in southern arabia 55,000 years ago. new evidence from the wadi surdud middle paleolithic site complex, western yemen[END_REF]) has been dated by optical stimulated luminescence (OSL) to 55 ka, providing insight into the Middle Paleolithic peopling of the Arabian Peninsula. The archeological layer is interstratified within thick sandy silt floodplain deposits filling a piedmont basin. Luminescence dates, lack of soil development, and gypsum precipitation indicate a high accretion rate of the floodplain during Marine Isotope Stage 3. The chronology of SD1 sequence (noted Yemen-sequence hereafter) was established by way of OSL dates determined on 498 single quartz grains (aliquots) from 16 sediment samples taken throughout the stratigraphy. Age determinations are based on Analyst software using the finite mixture model (FMM) [START_REF] Roberts | Distinguishing dose populations in sediment mixtures: a test of single-grain optical dating procedures using mixtures of laboratory-dosed quartz[END_REF][START_REF] Tribolo | Kobo 1 and l'abriaux-vaches (mali, west africa): two cases study for the optical dating of bioturbated sediments[END_REF]. Results for each sample are presented in terms of age and are obtained by dividing the equivalent aliquot paleodose by the dose rate (see Table 1 in supplementary material, p. 8, [START_REF] Delagnes | Inland human settlement in southern arabia 55,000 years ago. new evidence from the wadi surdud middle paleolithic site complex, western yemen[END_REF]. To convert the ages A in calendar dates, we apply a linear calibration curve such that t = t 0 -age, where the reference year t 0 is set equal to zero. It corresponds to the year of measurement at the laboratory: 2008/2009. The dates are positive and thus the prior time range T is set equal to [0 , 500] in order to largely include the data.

Following our modeling framework, each sample composed of a large set of dated quartz grains corresponds to one deposit Event (in other words all the aliquots in a sample refer to a same age) and the 16 Events are related to each other according to the stratigraphy, from the oldest sample Yemen-20b ( θ 1 ) to the youngest Yemen-2 (θ 16 ) (Fig. 14). Thus their true dates (Event) have to check the following relationship during the MCMC sampling:

θ 1 < θ 2 < ... < θ 16 .
The dispersion of the ages in each sample is important (see for instance Yemen-4 in Figures 17 and18) and could be explained by grains which were properly bleached at the time of deposit, or by some bioturbations which may occur and which are not easy to detect, or also by heterogenous alpha and beta dose rates in the sediment. This implies a lot of chrono-stratigraphic inversions as shown in Figure 15 when considering Events without any stratigraphic information. When the stratigraphic information is taken into account as in Figure 16, the Events are considerably improved and more precise. Extremal Events Yemen-20b and Yemen-2 are not changed but Events in between are shifted in order to respect the stratigraphy and HPD intervals become more precise as shown in Table 1. Typically, for Event Yemen-4, the 95% HPD age interval passes from [64.8, 75.5] without a stratigraphic relationship to [51.7, 54.6] with a stratigraphic order constraint.

Conclusion

The Bayesian Event model aims to estimate the date θ of a context (unit of stratification) from individual dates t i which are affected by errors of different types: laboratory and calibration curve errors and irreducible errors related to contaminations, taphonomic problems, etc, hence the possible presence of outliers. The Event approach is robust to outliers in the sense that individual variances σ 2 i put on all the dates t i in an Event act as outlier penalization. The Event model has a hierarchical structure which makes it possible to distinguish between date θ of an Event (in the sense of something that occurs in a certain place at a given time) and dates t i of the artifacts involved. Indeed, we assume that these artifacts, dated by different methods (chronometric methods, typo-chronology, history), are all contemporaneous to the Event. However some of them may be misallocated: the date can be reliable in itself because it corresponds to, for instance, the making of the artifact, but at the same time may not be contemporaneous to the target Event. The posterior distribution of the variance σ 2 i indicates if an observation is an outlier or not. However, it is not necessary to discard outliers because the corresponding high values σ 2 i will automatically penalize their contributions to the Event estimation. Our outlier model can apply without having to decide whether the outlier comes from a laboratory measurement error or from a dating error. Moreover, there are no exogenous or hyper parameters to adjust according to the different data of submitted sets. The only parameter involved in prior shrinkage, s 2 0 , comes uniquely from the data analysis via the calibration process. So, the approach is adapted to very different datasets. The good robustness properties of the Event model are paid with less precision in the dates. However, this loss of precision is compensated by better reliability of the chronology. Consequently, the Event model constitutes the basic element in our chronological modeling approach. In "Chronomodel" software, dating of a context, stratigraphic chronology or phasing are directly constructed by using the Event model.

Figure 5 .

 5 Figure5. Graphical representation of the densities defined in (15) and (17) with the following parameter values: s ji = 30 , σ gji = 10, σ δ = 10 2 and θ j = t j = 1000

Figure 6 .

 6 Figure 6. Synthetic data, Ex 1. [left] Event model: posterior densities t i (white background) and posterior density for Event θ (gray background). The individual posterior calibrated densities are superimposed in black. [right] t-type outlier model: posterior densities obtained with Oxcal software V4.2.

Figure 7 .

 7 Figure 7. Synthetic data, Ex 1. Event model: posterior densities of the individual standard deviations σ i in the Event model.

Figure 8 .

 8 Figure 8. Synthetic stratigraphic sequence with date inversion, Ex 2. Dates are Gaussian distributed with constant standard deviation equal to 30. Each Event contains 2 dates.

Figure 9 .

 9 Figure 9. Synthetic data in a stratigraphy with date inversion, Ex 2. [left] Event model: posterior densities t i (white background) and posterior density for Event θ (gray background). The individual posterior calibrated densities are superimposed in black. [right] t-type outlier model: posterior densities obtained with Oxcal software.

Figure 10 .

 10 Figure 10. Shroud of Turin, Ex. 3. [white background ] Posterior densities of t i and individual posterior calibrated densities (black line) obtained for 12 14 C dates. [gray background] Posterior density for Event θ.

Figure 11 .

 11 Figure 11. Turin, Ex. 3 (Cont.) Posterior densities obtained for standard deviations σ i .

Figure 12 .

 12 Figure 12. Lezoux Ex. 4. [white background ] Posterior densities of t i and individual posterior calibrated densities (black line) obtained for TL dates, for 14 C dates and for AM dates. [gray background ] Posterior density for Event θ

Figure 14 .

 14 Figure 14. Yemen sequence, Ex. 5. Stratigraphic model for the 16 Events, involving OSL dating on 498 single quartz grains.

Figure 15 .

 15 Figure15. Yemen sequence, Ex. 5. Posterior densities obtained for the events θ j , j = 1, ..., 16 when no stratigraphic information is put on the Events.

Figure 16 .Figure 17 .

 1617 Figure 16. Yemen sequence, Ex. 5. Posterior densities obtained for the events θ j , j = 1, ..., 16 with stratigraphic model defined in Fig. 14.

Figure 18 .

 18 Figure 18. Yemen, Ex. 5 (cont.). Posterior densities obtained for the standard deviations in Event Yemen-4.

This project is supported by the grant ANR-11-MONU-007 Chronomodel.

1