Lanos Philippe

Philippe Anne
email: anne.philippe@univ-nantes.fr.

EVENT DATE MODEL: A ROBUST BAYESIAN TOOL FOR CHRONOLOGY BUILDING

Keywords:

. We finally compare our prior on event dates implemented in ChronoModel with the prior in BCal and OxCal which involves supplementary parameters defined as boundaries to phases or sequences. Target event date model ; robust combination of dates ; prior information on dates ; hierarchical Bayesian statistics ; individual errors ; outlier penalization ;

Introduction

Bayesian chronological modeling appears as an important issue in archeology and palaeo-environmental sciences. This methodology has been developed since the 1990s [START_REF] Bayliss | Rolling out revolution: Using radiocarbon dating in archaeology[END_REF][START_REF] Bayliss | Quality in bayesian chronological models in archaeology[END_REF] and is now the method of choice especially for the interpretation of radiocarbon dates. 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[START_REF] Bronk Ramsey | Developments in the calibration and modelling of radiocarbon dates[END_REF][START_REF] Bronk Ramsey | Recent and planned developments of the program OxCal[END_REF].

All of these models provide an estimation of a chronology of dated events (DE), the estimated dates correspond to the dates of events that are actually dated by any chronometric technique. However calibrated radiocarbon dates and date estimates from other chronometric dating methods such as thermoluminescence, archaeomagnetism, dendrochronology, etc. can be combined with prior archeological information of various kinds to produce a combined chronology that should be more reliable than its individual components.

We propose a new chronological Bayesian model based on the concept of target event. This is related to the concepts of the dated event and the target event proposed by [START_REF] Dean | Independent dating in archaeological analysis[END_REF]. The target event (TE) is the event to which the date is to be applied by the chronometrician. Usually, the target events are not directly related to the dated events. According to [START_REF] Dean | Independent dating in archaeological analysis[END_REF], the dated event is contemporaneous with the target event when there is convergence (the two dates are coeval: DE = TE) and when the date is relevant to the target event. "Relevance" refers to the degree to which the date is applicable to the TE : it must be demonstrated or argued on the basis of archaeological or other evidence. Because it is generally not easy or possible to assure the relevance of the dated event to the target event date of interest, it is recommended to get many dates of dated event , and if possible from different dating techniques.

These dates can be outliers without having a means to determine if they are dating anomalies. In other words, we do not have any convincing archaeological arguments for rejecting them before modeling. This motivates the development of a robust statistical model for combining theses dates in such a way that it is very little sensitive to outliers (see [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. .

The target event model is a statistical model introduced in (Lanos and Philippe, 2017) for estimating the date of an event called "target event" . This model allows to combine in a robust way the dates of artifacts, which are assumed to be contemporary to this target event. To validate the robustness of this model to outliers, we provide a comparison of our model with the t-type outlier model implemented in Oxcal application. Numerical experiments also illustrate the sensibility to the outliers.

The target event model is then integrated to a most global model for constructing chronologies of target events. Prior information is brought upon the dates of the target events. We show how the prior archeological information based on relative dating between the target events (in a stratigraphic sequence for instance) or based on duration , hiatus or Terminus post quem (TPQ) or terminus ante quem (TAQ), assessment can be included in the model. The simulations and application illustrates the improvement brought in term of robustness. The proposed model is implemented in ChronoModel software, whose features are described in [START_REF] Vibet | Chronomodel v1.5 user's manual[END_REF]. This model can be compared with the standard chronological models when the target event is associated to only one dated event. However the interest of our approach is to get a robust estimation of the date of target event even if the dated events embeded dated events are outliers. We compare in the simulation part this approach with an outlier model based using discrete mixture distributions.

Many issues in archaeology raise the problem of phasing, that is how to characterize the beginning, the end and the duration of a given period. This question can be viewed as a post processing of the chronological model as for instance in Philippe and Vibet (2017a); [START_REF] Guérin | Chronology of the upper pleistocene loess sequence of havrincourt (france) and associated palaeolithic occupations: A bayesian approach from pedostratigraphy, osl, radiocarbon, tl and esr/u-series data[END_REF]. It is also possible to include additional parameters to characterise the phases. This requires the construction of a prior distribution on the dates of dated events, which belong to the phase. This point is discussed in this paper, we analyse the choice of model implemented in Oxcal application. We get an explicit form for the prior distribution of the dates, and we show that this probability distribution behaves in the same way as the dates included in a target event in the sense that we observe a concentration of the dates belonging to the phase. Our result also shows the difficulty to construct appropriate prior on a sequence of dates.

The paper is organised as follows : In Section 2 we recall the construction of the target event model and we provide a theoretical comparison with the t-outlier model. In Section 3 we propose a new Bayesian model to estimate dates of target event. In Section 4, we analyse the Bayesian modelling of phases implemented in Oxcal application, and we compare this model with the target event model. In Sections 5 -6 we apply our models to simulated and real datasets.

A robust combination of dates : the target event date model

We propose a statistical Bayesian approach for combining dates, which is based on a robust statistical model. This combination of dates is aiming to date a called "target" event defined on the basis of archaeological/historical arguments.

2.1. Description of the Bayesian model. We propose to use a hierarchical Bayesian model to estimate the date θ of a target event Et. It combines dates t i (i = 1, ..., n) of dated events Ed in such a way that it is robust to outliers. This model is based on very few assumptions and does not need to tune hyperparameters. It is described as follows.

The event dates t i are estimated from n independent measurements (observations, also called determinations in Buck's terminology) M i yielded by the different chronometric techniques. Each measurement M i , obtained from a specific dating technique, can be 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.

At this step, the random effect model can be written as follows:

M i = µ i + s i i , ∀ i = 1, ..., n µ i = g i (t i ) + σ gi (t i )ρ i (1)
where ( 1 , ... n , ρ 1 , ..., ρ n ) are independent and identically Gaussian distributed random variables with zero mean and variance 1, and where :

• s i i represents the experimental error provided by the laboratory.

• σ gi (t i )ρ i represents the error provided by the calibration curve. We assume that the experimental error provided by the laboratory are independent. Dependence structure could be added between the (s i i ) i=1,...,n in order to take into account the systematic error associated to the dating techniques of each laboratory (see [START_REF] Combs | Bayesian analysis of individual and systematic multiplicative errors for estimating ages with stratigraphic constraints in optically stimulated luminescence dating[END_REF] in the particular case of the optically stimulated luminescence dating). However the construction of such a model requires additional information depending on each dating method and each laboratory, which is not available in most of the application of chronological models.

Remark 1. 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 according to the R-combine model [START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF] before being incorporated in the event date model. This model can be viewed as a degenerated version of the target event model by taking σ i = 0 for all i. On the contrary, if several calibration curves are involved, the R-combine model is no longer valid and chronometric dates are directly incorporated in the event date model.

The target event date θ is estimated from the dated events t i . The main assumption we make is the contemporaneity (relevance and convergence, following [START_REF] Dean | Independent dating in archaeological analysis[END_REF]), of the dates t i with the event date θ. However, because of error sources of unknown origin, it can exist an "over-dispersion" or dating anomalies of the dates with respect to θ. These errors can come from the way of ensuring that the samples studied can realistically provide results for the events that we wish to characterize, the care in sampling in the field, the care in sample handling, preparation and measure in the laboratory, or other non-controllable random factors that can appear during the process [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF]. This over-dispersion is also similar to the irreducible error described in [START_REF] Niu | The Bayesian approach to radiocarbon calibration curve estimation: the intcal13, marine 13, and shcal13 methodologies[END_REF] in the framework of radiocarbon calibration curve building.

Consequently, we model the over-dispersion by an individual error σ i according to the following random effect model:

(2)

t i = θ + σ i λ i
where (λ 1 , ..., λ n ) are independent and identically Gaussian distributed random variables with zero mean and variance 1. The individual error σ i measures the degree of disagreement which can exist between a date (Ed) and its target event date (Et). It will be a posteriori small when the dates (Ed) are consistent with the target date θ and high when the date t i is far from the target date θ. Consequently, the posterior distribution of the individual error σ i will give some information about the "outlying" state of a date with respect to the target date. Finally, the joint distribution of the probabilistic model can be written according to a Bayesian hierarchical structure: (3)

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 ), σ 2 i ∼ Shrink(s 2 0 ). (4)
The parameter of interest θ is assumed to be uniformly distributed on an interval

T = [T a , T b ] (5) θ ∼ Uniform(T ).
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 imparted to the dates t i . Consequently their support is the set of real numbers R.

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

(6) p(σ 2 i ) = s 2 0 (s 2 0 + σ 2 i ) 2 1 [0,∞] (σ 2 i ),
where 1 A (x) is the indicator function (= 1 if x ∈ A, = 0 if x / ∈ A) and where the parameter s 2 0 must be fixed. The motivation for this choice of prior is described in detail in [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF]. Parameter s 2 0 quantifies the magnitude of error on the measurements. It is estimated according to the following process:

(1) An individual calibration step is done for each measurement M i , i = 1, ..., n.

It consists of the simple model

M i |t i ∼ N (g i (t i ), s 2 i + σ 2 gi (t i )), t i ∼ Uniform(T ).
using the same notation as (4) and (5). For each i, (a) Sample from the posterior distribution of

t i given M i (7) p(t i |M i ) ∝ 1 S i exp -1 2S 2 i (M i -g i (t i )) 2 1 T (t i ),
where

S 2 i = s 2 i + σ 2 gi (t i ) (b) Approximate the posterior variance var(t i |M i ) by its Monte Carlo ap- proximation denoted w 2 i .
(2) Take as shrinkage parameter s 2 0 : 1

s 2 0 = 1 n n i=1 1 w 2 i .
Remark 2. The calibration curve g i in radiocarbon or in archaeomagnetic dating is always defined on a bounded support [T m , T M ]. Therefore an extension of g i on R is required to well defined the conditional distribution of t i in (4).

We suggest to extend the calibration curve by an arbitrary constant value with a very large variance with respect to the known reference curve: for instance

g i (t) = (g i (T m ) + g i (T M ))/2 and σ 2 gi = 10 6 sup t∈[Tm;T M ] (g i (t)) - inf t∈[Tm;T M ] (g i (t)) 2
In the case of TL/OSL or Gauss measurements, there is no need for an extension because g i (t) is defined whatever t in R.

This statistical approach does not model the outliers, i.e. we do not estimate the posterior probability that a date is an outlier. Outlier modeling (See section 2.2 for such an approach) can provide more accurate results, but it often requires two (maybe more) estimations of the model: 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 (however, the question remains to know when to stop, because new outliers can appear during the subsequent runs). The event model which is based on the choice of robustness, avoids this two-step procedure.

2.2.

Comparison with alternative outlier event models. Several outlier models have been implemented since the 90's, in the framework of radiocarbon dating. So we can distinguish between mainly three outlier modelings:

(1) Outlier model with respect to the measurement parameter [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF][START_REF] Buck | Bayesian tools for tephrochronology[END_REF].

(2) Outlier model with respect to the laboratory variance parameter [START_REF] Pérez | A new robust statistical model for radiocarbon data[END_REF]. ( 3) Outlier model with respect to the time parameter. This corresponds to the t-type outlier model [START_REF] Bronk Ramsey | Dealing with outliers and offsets in radiocarbon dating[END_REF] The t-type outlier model can be compared to the event model, by considering that the date parameter t directly coincides with our target event date parameter θ. It means that the hierarchical level between t and θ does no longer exist in the t-type outlier model. Hence the observed measurement will be noted M j instead of M i , the true unknown measurement noted µ j instead of µ i and a date noted t j instead of t i in order to avoid any confusion with index i used for chronometric dates (dated events) in the target event date model.

The model with random effect becomes :

M j = µ j + s j j ,
µ j = g j (t j + δ j φ j 10 u ) + σ gj (t j + δ j φ j 10 u )ρ j (8) where

• ( 1 , ... r , ρ 1 , ..., ρ r ) are independent and identically Gaussian distributed random variables with zero mean and variance 1. • s j j represents the experimental error provided by the laboratory and σ gj (t j + δ j φ j 10 u )ρ j the calibration error. • the prior on φ j is the Bernoulli distribution with parameter p j . A priori, φ j 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 corresponds to the shift on the measurement M j if it is detected as an outlier. • u is a scale parameter to offset δ j . Parameter u can be fixed (for instance 0) or prior distributed as Uniform([0,4])

When parameter u is fixed, the joint distribution of the probabilistic model for r dates t j (j = 1, ..., r) can be written in the form:

(9) p(M, µ, t, δ, φ) = r j=1 p(M j |µ j )p(µ j |t j , δ j , φ j )p(t j )p(δ j )p(φ j ),
where the conditional distributions that appear in the decomposition are given by:

M j |µ j ∼ N (µ j , s 2 j ) µ j |t j ∼ N (g j (t j + δ j φ j ), σ 2 gj (t j + δ j φ j )) t j ∼ Uniform([T a , T b ]) δ j ∼ N (0, σ 2 δ ) or ∼ T (ν) : Student's t-distribution with ν degrees of freedom φ j ∼ Bernoulli(p j )
This modeling can be compared to the event date model if we consider several dates nested in the Oxcal function "Combine" with t-type outlier model. The dates t j (j = 1, ..., r) then becomes a common date t. Consequently, ( 9) is transformed into:

(10) p(M, µ, t, δ, φ) ∝ p(t) r j=1 p(M j |µ j )p(µ j |t, δ j , φ j )p(δ j )p(φ j ),
If we consider linear calibration curves with constant errors, that is, by setting:

• g j (t) = t • σ gj (t) = σ g ,
and knowing that

• δ j ∼ N (0, σ 2 δ ) • φ j ∼ Bernoulli(p),
it is possible to analytically integrate (10) with respect to δ j , φ j , and µ j . The posterior probability density of t is then given by:

(11) p(t|M ) ∝ p(t) r j=1 p(M j |t),
where the conditional distribution of M j given t is a finite mixture distribution defined by ( 12)

p(M j |t) = p 1 √ 2π s 2 j + σ 2 g + σ 2 δ e -(M j -t) 2 2(s 2 j +σ 2 g +σ 2 δ ) + (1 -p) 1 s 2 j + σ 2 g √ 2π e -(M j -t) 2 -2(s 2 j +σ 2 g ) .
Note that the outlier model operates well as a model averaging thanks to the structure of mixture distribution in (12).

On the other hand, posterior density of event date θ = t deduced from ( 3) can be compared to density of time t in (11) after integration with respect to µ i = µ j and σ 2 i = σ 2 j . The posterior probability density of θ is given by:

(13) p(θ|M ) ∝ p(θ) r j=1 p(M j |θ),
where ( 14)

p(M j |θ) = ∞ 0 1 2π(s 2 j + σ 2 g + σ 2 j ) e 1 -2(s 2 j +σ 2 g +σ 2 j ) (Mj -θ) 2 s 2 0j (s 2 0j + σ 2 j ) 2 dσ 2 j
A graphical representation of the densities defined in ( 12) and ( 14) are shown in Figure 1 with the following parameter values: s j = 30 , σ g = 10, σ δ = 10 2 and θ = t = 1000. The density ( 12) is plotted for three different values p = 0.01, 0.05, 0.10 . We can observe that shrinkage modeling in ( 14) leads to a more diffuse density making it possible to better take into account the possible presence of outliers. 12) and ( 14) with the following parameter values: s j = 30 , σ g = 10, σ δ = 10 2 and θ = t = 1000

Prior information on groups of target event dates

We now consider a group of dates θ j (j = 1, ..., r) of target events which belong either to a "stratigraphic phase defined as a group of ordered contexts, or to a "chronological phase defined as a set of contexts build on the basis of archaeological, architectural, geological, environmental,...criteria. In practice, a "context" is defined by the nature of the stratification at a site, and the excavation approach used by the archaeologist. Together, these two determine the smallest units of space and time -the context-that can be identified in the stratigraphic record at an archaeological site. Because the term "phase" is also used to describe Bayesian chronological models, [START_REF] Dye | Archaeological sequence diagrams and bayesian chronological models[END_REF] use the term "stratigraphic phase to refer to a group of ordered contexts, and the term "chronological phase to refer to a time period in a chronological model.

The concept of phase is currently implemented in chronological modeling by using a specific parametrization which will be discussed in details in Section 4. Hereafter, we extend the target event date model in a simple way to the case of a group of dates related to each other by order or duration relationships. Hence we propose to estimate the beginning, the end and the duration of a phase directly from the group of target dates, without adding any supplementary parametrization.

3.1. Characterization of a group of target event dates. We consider a collection of r event dates denoted θ j with j = 1, ..., r. For each event date θ j , we suppose that n j measurements denoted M j = (M j1 , . . . , M jnj ) charactering the event are observed. We assume the independence of the measurements conditionnally to the event dates (θ 1 , ..., θ r )

p(M 1 , . . . , M r , θ 1 , ..., θ r ) = p(θ 1 , ..., θ r ) r j=1 p(M j |θ j )
where p(M j |θ j ) is the conditionnal distribution of M j given θ j in the hierarchical model defined in [START_REF] Vibet | Bayesian Modeling of Archaeological Chronologies[END_REF]. The Bayesian model is summarized by Figure 2.

When there is no supplementary information, we assumed that the dates θ j are independent and uniformly distributed on the time interval (study period) T = [T a , T b ]. We remind that this interval is fixed by the user based on historical or archeological evidences. Consequently, the prior on vector θ is such that:

(15) p(θ 1 , ..., θ r ) = 1 (T b -T a ) r r j=1 1 T (θ j )
The overall Bayesian model is described in Figure 2 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. 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). Almost always we do not know how a group of event dates can a priori be distributed in a phase. It means that, in our approach, a phase does not respond to a statistical model. From a chronological point of view, all the information is carried by the target event dates themselves. Consequently, we estimate the beginning and the end of a phase as a measurable application of the parameters θ j (j = 1, ..., r). Thus the beginning of a phase is estimated by the minimum θ (1) of the r event dates included in the phase ( 16) θ (1) = min(θ j , j = 1, . . . , r).

j = 1 to r i = 1 to n j θ j t ji µ ji M ji σ 2 ji s 2 ji
In the same way, the end of a phase is estimated by the maximum θ (r) of the r event dates included in the phase: (17) θ (r) = max(θ j , j = 1, . . . , r).

By the plug-in principle, we estimate the duration of the phase by

(18) τ = θ (r) -θ (1)
Considering two phases P k = {θ k,1 , ..., θ k,r k }, k = 1, 2, the hiatus between P 1 and P 2 is the time gap between the end of P 1 and the beginning P 2 . The hiatus is estimated by :

(19) γ = max(θ 2,(1) -θ 1,(r1) , 0)
where θ 2,(1) = min(θ 2,1 , ..., θ 2,r2 ) and where θ 1,(r1) = max(θ 1,1 , ..., θ 1,r1 ). Note that the estimate γ takes the value 0 if θ 1,(r1) ≤ θ 2,(1) , that corresponds to the absence of hiatus. The conditional distributions of the parameters θ (1) , θ (r) , γ and τ given the observations can be easily derived from the joint posterior distributions of the event dates.

Remark 3. Let us remind that these parameters are estimated knowing the data, that is to say from the target event dates which are available in the phase of interest.

In particular, a good estimation of the beginning (resp. the end) of a phase requires that the archaeologist has sampled artifacts belonging to target events which are very near to this beginning (resp. this end).

These estimates are valid whatever the prior on the event dates are. Precision on the estimation of these parameters can be gained if it is possible to add supplementary information on the target event dates within the study period T . For instance some temporal orders induce some restrictions on the support of the distribution.

Different prior on dates θ j can be defined according to the following circumstances:

• 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], can imply antero-posteriority relationships between target dates θ j . This can also imply antero-posteriority relationships between groups of dates θ j , these groups defining different phases. • We can have some prior information about the maximal duration of a group of target events. • We can also have some prior information about the minimal temporal hiatus between two groups of target events. It has been demonstrated that these types of supplementary prior information can significantly improve chronometric dates (see initial works of [START_REF] Naylor | An archaeological inference problem[END_REF]; [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 the next three subsections, we discuss each of these prior information in the framework of the target event date model. Our modeling approach is very soft thanks to the fact that all the relationships operate directly onto the target event dates θ j . As a consequence, in a global modeling project, different phasing systems (multiphasing) can be defined on the basis of different criteria and so a phasing system can intersect an other phasing system: for example a ceramic phasing can intersect a lithic phasing in the sense that some target events can belong to two or several phases. Such a constraints network, when available, can significantly contribute to improve the precision of the estimates.

3.2. Prior information on temporal order. Target events or phases (groups) of target events can have to check order relationships. This order can be defined in different ways: by the stratigraphic relationship (physical relationship observed in the field) or by stylistic, technical, architectural,... criteria which may be a priori known. Thus, the constraint of succession is equivalent to a hiatus of unknown amplitude put between event dates or groups (phases) of event dates.

If we consider a stratigraphic sequence composed of target events, the prior on vector θ becomes:

(20) p(θ 1 , ..., θ r ) ∝ 1 (T b -T a ) r 1 C (θ 1 , ..., θ r ) with C = S ∩ T r and • T r = [T a , T b ]
r the support which defines the study period,

• S = the group of r-uplets event dates θ j which respect total or partial order relationships. We can also consider two groups of target events, P k = {θ k,1 , ..., θ k,r k }, k = 1, 2, containing r k event dates such that all the dates of P 1 are before all the dates of P 2 . The following equations should be verified by all the events included in both phases.

∀j ∈ {1, ..., r 1 }, ∀l ∈ {1, ..., r 2 }, θ 1,j < θ 2,l or max(θ 1,j , j = 1, . . . , r 1 ) < min(θ 2,l , l = 1, . . . , r 2 ) or ( 21) min(θ 2,l , l = 1, . . . , r 2 ) -max(θ 1,j , j = 1, . . . , r 1 ) > 0

The r-uplets event dates θ j will then have to check total or partial stratigraphic constraints and also to satisfy the inequality (21) between event dates of the two phases. We can see that a succession constraint between phases operates in the same mathematical way as for a set of stratigraphic constraints put between all the individual target events.

Remark 4. Estimating the date of a target event needs to incorporate several dates (Ed), otherwise the event date modeling will not yield a better posterior information. However, it is possible to nest only one date per target event provided that the group of events is constrained by temporal order.

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 target event date 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. This is illustrated in Section 5.2.

3.3.

Prior information about the duration. Prior information can be included on the duration of a phase, that is on a group of target events. We can impose a maximal duration τ 0 . This means that all the event dates θ j (j = 1, ..., r) in the phase have to verify the constraint of duration according to the following equation: (22) max(θ j , j = 1, . . . , r) -min(θ j , j = 1, . . . , r) ≤ τ 0

The r-uplets event dates θ j will have to check total or partial stratigraphic constraints and also to satisfy the inequality 22 between event dates in the phase. This means that any r-uplet θ j sampled during the MCMC process has a duration of at the most τ 0 .

3.4. Prior information about the amplitude of a hiatus. Prior information about a hiatus γ between two phases P k = {θ k,1 , ..., θ k,r k }, k = 1, 2, may be available. Hence we can impose that the amplitude of the hiatus is higher than a known value γ 0 > 0. All the event dates of phase P 1 and phase P 2 have to verify the following constraint :

(23) min(θ 2,l , l = 1, ..., r 2 ) -max(θ 1j , j = 1, ..., r 1 ) ≥ γ 0

The event dates will have to check total or partial stratigraphic constraints and also to satisfy the inequality ( 23). This means that any r-uplet θ 1,j from phase P 1 and sampled during the MCMC process is separated by a time span of at least γ 0 from the r-uplet θ 2,l of the next phase. Note that it is obviously not possible to impose a hiatus between two phases when a same event belongs to these two phases.

3.5. Prior information about known event dates: the bounds. Bounds, such as historical dates, Terminus post quem (TPQ) or terminus ante quem (TAQ), may also be introduced in order to constrain one or several event dates θ j . If we consider a set of r events assumed to happen after a special event with true calendar date B such that B < (θ 1 , ..., θ r )

This condition must be included in the set of constraints that define the support of the prior distribution of the event dates. The bound can be also defined with an uncertainty, i.e. B ∈ [B a , B b ] ⊂ [T a , T b ], and so it is included in the set of the parameters. The prior density of the event dates can then be written as:

p(θ 1 , ..., θ r , B) = p(θ 1 , ..., θ r |B)p(B)
where

p(θ 1 , ..., θ r |B) = 1 (T b -B) r r j=1 1 [B,T b ] (θ j ) and p(B) = 1 B b -B a 1 [Ba,B b ] (B)
It is important to note that the introduction of bounds B in a global model composed of groups (phases) of event dates must remain consistent with other constraints of temporal order, of duration or of hiatus described in previous sections.

Discussion on date prior probabilities

In this section we discuss the question of the specification of the prior on the event dates θ j . When there is no supplementary information, it may seem natural to assume that the dates θ j are independent and uniformly distributed on the time interval (study period) T = [T a , T b ]. However, from an archaeological point of view, it may seem more natural to assert that the span ∆ = max(θ j ) -min(θ j ) has a uniform prior distribution We can see that this assumption is not checked when starting from the prior density (15). After an appropriate change of variable and letting R = T b -T a , we can determine the density of the span according to the number of dates r and regardless their order:

(24) p(∆) = r(r -1) R r (R -∆)∆ r-2
A span of 2∆ is favoured over a span of ∆ by a factor of about ∆ r-2 when ∆ R (spreading tendency when r becomes large). This behavior occurs regardless of the order between dates. [START_REF] Naylor | An archaeological inference problem[END_REF]; [START_REF] Buck | Calibration of radiocarbon results pertaining to related archaeological events[END_REF]; [START_REF] Christen | Summarizing a set of radiocarbon determinations: a robust approach[END_REF]; [START_REF] Buck | The Bayesian Approach to Interpreting Archaeological Data[END_REF]; [START_REF] Nicholls | New radiocarbon calibration software[END_REF] propose a specific phase modeling which aims to avoid this spreading bias. We study the properties of this modeling implemented in BCal and OxCal software. The Naylor-Smith-Buck-Christen (NSBC) prior is defined for a group of event dates which are placed between two additional hyperparameters α and β in the Bayesian hierarchical structure. These hyperparameters α and β represent boundaries where α is the beginning and β the end of the group of events called Phase. The dates θ j are assumed to be independent conditionally to these two boundaries. Then, the prior density on vector θ becomes:

(25) p(θ 1 , . . . , θ r |α, β) = 1 (β -α) r 1 [α,β] r (θ 1 , . . . , θ r )
In the absence of supplementary information, a non informative prior density is assigned to α and β : From these hypotheses, it is possible to calculate the prior joint probability density of the event dates (θ 1 . . . θ r ). This is carried out by integration against α and β between the two limits T a and T b . For 2 (ordered or non ordered) event dates θ 1 , θ 2 , we obtain:

p(θ 1 , θ 2 ) = p(θ 1 , θ 2 |α, β).p(α, β) dαdβ = 1 (β -α) 2 .1 [α,β] (θ 1 ).1 [α,β] (θ 2 ) 2 (T b -T a ) 2 .1 [Ta≤α≤β≤T b ] (α, β). dαdβ = 2 (T b -T a ) 2 min(θ1,θ2) Ta T b max(θ1,θ2) 1 (β -α) 2 dβ dα = 2 (T b -T a ) 2 min(θ1,θ2) Ta [ 1 (max(θ 1 , θ 2 ) -α) - 1 (T b -α) ] dα
Finally, we have:

(27) p(θ 1 , θ 2 ) = 4 (T b -T a ) 2 (-ln(max(θ 1 , θ 2 ) -min(θ 1 , θ 2 )) + ln(T b -min(θ 1 , θ 2 )) + ln(max(θ 1 , θ 2 ) -T a ) -ln(T b -T a ))1 [Ta,T b ] 2 (θ 1 , θ 2 )
The joint probability density for r (ordered or non ordered) event dates (θ j , j = 1, . . . , r) with r ≥ 3, is obtained by following the same integration process:

(28) p(θ 1 , ..., θ r ) = 2 r! (r -1)(r -2)(T b -T a ) 2 1 [θ (r) -θ (1) ] r-2 - 1 [T b -θ (1) ] r-2 - 1 [θ (r) -T a ] r-2 + 1 [T b -T a ] r-2 1 [Ta,T b ] 2 (θ (1) , θ (r) )
where θ (r) = max(θ j , j = 1, . . . , r) and θ (1) = min(θ j , j = 1, . . . , r).

Equations ( 27) and ( 28) clearly show that the priors ( 25) and ( 26) provoke a strong concentration effect of the dates θ j as illustrated in Figure 3 which is calculated for two event dates with formula ( 27). The region of high probability of (θ 1 , θ 2 ) is concentrated around the first diagonal (θ 1 = θ 2 ). In conclusion, the NSBC prior clearly favors the fact that dates (θ 1 , ..., θ r ) become near each other. This property of the phase can be compared with the assumption of contemporaneity imposed on the dates in the event model. Starting from the prior densities ( 27) and ( 28), and after an appropriate change of variable, we can also determine the density of the span ∆ = max(θ j , j = 1, ..., r)min(θ j , j = 1, ..., r), according to the number of dates r. Letting R = T b -T a , we obtain :

r = 2 : p(∆) = 4 R 2 (R + ∆)(ln(R) -ln(∆)) + 2(∆ -R)) (29) r = 3 : p(∆) = 6 R 2 (R -∆)(1 + ∆ R ) + 2∆(ln(∆) -ln(R)) (30) r ≥ 4 : p(∆) = 2r (r -2)R 2 (R -∆)(1 + ∆ r-2 R r-2 ) - 2 r -2) (∆ - ∆ r-2 R r-3 ) (31)
These distributions are shown in Figure 4. When r is small (r = 2, 3, 4), the densities are very high for ∆ near zero, hence a high concentration effect. When r tends towards infinity, the distribution of ∆ tends to p(∆) = 2 R 2 (R -∆). This latter triangular distribution is still maximal when ∆ = 0, and it becomes more and more flat as R increases. [START_REF] Nicholls | New radiocarbon calibration software[END_REF] proposed to make the prior distribution on the dates θ j uniform by multiplying the prior in (26

) by R 2 2(R-∆)
where min(θ j , j = 1, ..., r) is replaced by α and max(θ j , j = 1, ..., r) by β . This option is implemented in OxCal when setting UniformSpanPrior=true' in a modeling project (Bronk Ramsey, 2009a). Instead of looking at the span ∆, we can look at the variance of the event dates θ j , which is more representative of their scattering. This variance is proportional to the sequence of the Euclidian distance of the event dates to the straight line θ 1 = θ 2 = ... = θ r in the space of dimension r and therefore it gives a good way to characterize the scattering of the dates. We observe that this distance remains near zero whatever the number r of event dates. Figure 5 shows an evaluation of the dispersion of the dates through the density of the statistic

(32) 1 r r i=1 θ i - 1 r r j=1 θ j 2
obtained with r = 10. Conditionally to (α, β), the event dates are uniformly sampled between α and β and the distribution of (α, β) is given by

p(α) = 2(T b -α) (T b -T a ) 2 p(β|α), = 1 (T b -α) 1 [α,T b ] (β),
Note that we simulate a sample from the prior distribution of (α, β) as follows :

1. Generate u ∼ Uniform(0, 1) 2. Take α = T b -R √ 1 -u 3. Generate β ∼ Uniform(α, T b ).
The same result is obtained whatever the partial or total order of the event dates. The concentration effect is particularly visible when considering dates with different uncertainties. In this case, the posterior results for α and β are attracted by the most precise dates.

Our result shows that the prior on event dates θ j given by formulas ( 27) and ( 28) does not give an ideal solution for the specification of the prior density on dates θ j . It allows the correction of the date spreading bias when the uncertainties on the dates are of the same amplitude, but it can generate some undesirable concentration effect when these uncertainties are different from each other. This behavior is illustrated in the following example.

Example 1. Synthetic data in NSBC phase

We consider a phase with 5 non ordered Gaussian dates: 3 dates 0 ± 30 AD, 1 date -500 ± 200 AD and 1 date +500 ± 200 AD. We can see in Figure 6 [left] that the posterior densities for α and β calculated with OxCal program are very close to the precise dates 0 ± 30, a result predicted by formula (27) . It is important to note that, in this case, the agreement index A given by OxCal are low and produce a warning to prevent from an over-interpretation of the result. Figure 3 [right] shows the result obtained with ChronoModel program when the five dates are nested in the event date model. We observe that the NSBC result tends to behave like the event date model in the sense that it favors the contemporaneity of the dates. In conclusion, the specification of an appropriate prior on target event dates still remains an open question. We showed that the NSBC phase model favors, in an underlying way, the concentration of the event dates in a phase. We also showed that the uniform prior favors spreading dates, especially when temporal order constraints act.

Simulations

In this section, Monte Carlo experiments are done to illustrate the performance of the proposed model on simulated samples. All the simulations are done using the R package ArchaeoChron (see Philippe and Vibet (2017c)).

5.1. The target event model. To illustrate the properties of the event date modeling, especially its robustness, we simulate sample of measurements M i with outliers using a mixture with two Gaussian components:

(33)

(1 -q)N (0, 1) + qN (µ, 1)

The date of the target event is 0, and q represents the proportion of outliers. We compare the event target model with • the r-combine model. See Remark 1.

• the t-outlier model described in Section 2.2 We assume that all the measurements have the same calibation curve, this condition is required for the r-combine model. We take g(x) = x and σ g = 0. The time range (study period) T is set equal to [-20 , 20], Figure 7 and 8 represent the boxplot of the Bayes estimates evaluated on 500 independent replications. This confirms that the r-combine model is not robust to outliers. For q = 0 (no outlier) the three models gives the same results, we do not observe a significant difference in terms of accuracy. The target event model ant the t-outlier model behave similarly, and provide robust statistical methods. Indeed for q < 5% the presence of outliers has practically no influence on the target date result in the sense that the boxplots stay centered around the true value

The lengths of the boxplots also indicate that the loss of accuracy if less important for event model than the t-outlier model. To understand the robustness of the target event model, we represent in Figure 8 the boxplot of the Bayes estimate of the individual variance σ 2 i . We compare the values of this parameter when the date t i of dated event is an outlier or not. Recall that the date t i is an outlier when it is not contemporaneous of the date θ. Figure 8 shows that the individual standard deviation takes large values for the outliers. This penalizes the contribution of the observation on the estimation of the date, leading to robust procedure. Moreover regardless to the proportion of outliers, the behavior is the same. The parameters are the same as in Figure 7 5.2. Chronological model. In this part, we illustrate the behavior of our chronological model in presence of stratigraphic inversion. We want to date target events, which satistfy stratigraphy constraints θ 1 < θ 2 < • • • < θ 10 . To compare our model with the t-outlier model, we construct target event containing only one dated event.

We denote t 1 , t 2 , . . . , t 10 the dates of dated event. To evaluate the robustness of event model, we assume that the date t 5 creates a stratigraphic inversion. In other words we take θ i = t i for all i = 5 and t 5 > θ 5

We simulate sample of 10 measurements such that the true dates satisfy the condition t 1 < t 2 < • • • < t 10 . Then, we shift the observation associated with the date t 5 to create the stratigraphic inversion in the observations. The prior information on the dates of target events is θ 1 < • • • < θ n while this constraints is imposed on the dates t 1 < t 2 < • • • < t 10 for the t-outlier model. The time range (study period) T is set equal to [-10 , 100]. Figure 9 and 10 compares our model with the t-outlier model.

Figure 9 gives the results for the target event model. The value of t 5 clearly appears as an outlier, as shown by the high values of the standard deviation σ 5 (see Fig 9). The estimations of t 5 are centered around the true value since the stratigraphic constraints is not imposed to the parameters t i . Due to the stratigraphic constraint, the Bayes estimate of θ 5 is far from the true values of t 5 . But this outlier does not disturb the estimation of the rest of the sequence. For all the other dates the stratigraphic constraints on θ i bring prior information and make more precise the target dates θ i than the dates t i .

Figure 10 gives the results of the t-type outlier model is applied to the same datasets for different values of the parameter σ δ . When σ δ = 0, the outlier are not taken into account in the model. For a good choice of the parameter σ δ (here 10), the performance of the t-outlier model is quite similar to those of the target event model. But Figure 10 shows that the t-type outlier model is very sensitive to the choice of its hyperparameter. The lack of an adaptative choice for this parameter is the main drawback of this approach. q q q q q q q q q q q q q q theta[1] theta[5] theta[9] 0 5 10 15 20 25 target event q q q q q q q q q q q q q q q q q q q q q q q t[1] t [START_REF] Vibet | Bayesian Modeling of Archaeological Chronologies[END_REF] 0 5 10 15 20 25 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q sigma[1] sigma[5] sigma q q q q q q q q q theta[1] theta [7] 0 5 10 15 20 25 q q q q q q q q q q q q q q q q q q q q dated event whitout outlier model q q q q q q q q q theta[1] theta [7] 0 5 10 15 20 25 q q q q q q q q q q q q q q q q q q q q dated event with outlier model sigma_delta = 1 q q q q q q q q q q q q q q q theta[1] theta [7] 0 5 10 15 20 25 q q q q q q q q q q q q q q q q q q q q dated event with outlier model sigma_delta = 10 q q q q theta[1] theta [7] 0 5 10 15 20 25 q q q q q q q q q q q q q q q q q q q q dated event with outlier model sigma_delta = 100 Figure 10. Boxplot of the Bayes estimates of t i for the t-outlier model with different values of the hyperparameter in (12) : σ δ varies and p = 5%. The simulated sample are the same as in Fig. 9 The circles (resp. stars ) indicate the true values of the dates θ i . (resp t i )

Applications

The model is implemented in the cross-platform application Chronomodel which is free and open source software [START_REF] Lanos | Chronomodel : Chronological Modelling of Archaeological Data using Bayesian Statistics[END_REF]. Features of the application are described in details in [START_REF] Vibet | Chronomodel v1.5 user's manual[END_REF] : graphical user interface for importing data and modeling construction, MCMC options and controls, graphical and numerical results. Different graphical tools are implemented for assessing the convergence of the MCMC: the history plot, the autocorrelation function, the acceptance rate of Metropolis-Hastings algorithms (see Appendix, for details) The user can adjust the length of burn-in, the maximum number of iterations for adaptation and for acquisition, and the thinning rate [START_REF] Vibet | Chronomodel v1.5 user's manual[END_REF].

6.1. Lezoux (Auvergne, France): last firing of a potter's kiln. The aim is to date the last firing of a medieval potter's kiln recovered at the 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). The target event corresponds to the last firing of the kiln, and this corresponds to the last use of the kiln (i.e. behavioral event, according to [START_REF] Dean | Independent dating in archaeological analysis[END_REF]). Some events relevant to this target event are dated by three chronometric techniques: archaeomagnetism (AM), thermoluminescence (TL) and radiocarbon (14C). AM and TL dates are determined from baked clay and 14C date from charcoals of trees assumed to be felled at the same time as the last firing.

The observations are composed of 3 TL dates (CLER 202a,202b,203), 2 AM dates (inclination and declination) and 1 radiocarbon date (Ly-5212) As the kiln belongs to the historical period, the prior time range T is set equal to [-1000 , 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. The target date gives a 95% HPD interval equal to [574,885] AD. Some of the posterior densities for standard deviations σ i (Fig. 12) are spread. 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 some favorable archaeological contexts, it is possible to get information about the duration of a phase. It is the case here with an household cluster excavated on the Malpaís Prieto site (Michoacán State, Mexico), in the framework of the archaeological project Uacusecha [START_REF] Pereira | Ephemeral cities? the longevity of the postclassic tarascan urban sites of zacapu malpas and its consequences on the migration process[END_REF]. The archaeological artifacts are typical of a chrono-ceramic phase which is dated to the period 1200 -1450 AD, so the prior time range (study period) T can be set equal to this interval. Five radiocarbon ages have been obtained from burials (bone samples) and a midden (charcoal samples). Looking at the individual calibrated radiocarbon dates, which appear to be consistent between them (no outliers), the overall date range for the occupation is between 1276 and 1443 AD (at 95% confidence level). Each radiocarbon date corresponds to one target event and these events allow to estimate the beginning and the end of the household cluster phase to which they belong to (Fig. 13,left). This modeling with only one date per target event and without any stratigraphic or duration constraints gives large 95% HPD intervals for the beginning and end estimates: Thus, this simple modeling does not allow to significantly improve the prior archaeological information. To do this, various archaeological evidences can better constrain the occupation duration: stratigraphic evidences, accumulation processes of the occupation remains, durability of the partly perishable architecture. In this example, these information indicate that the occupation has hardly exceeded one century, the most plausible estimation being between 60 and 90 years. As a trial, we consider that the phase duration cannot exceed τ 0 = 60 years. Applying this constraint during the MCMC sampling according to (22), we obtain the following results (Fig. 13 Posterior densities of the target event dates θ j in Figure 13 [right] are greatly shrunken compared with the unconstrained modeling. Moreover, this modeling tends to favor the recent part of the study period, by mitigating the bimodal shape of the individual calibrated dates.

Conclusion

The Bayesian event date model aims to estimate the date θ of a target event from the combination of individual dates t i coming from relevant dated events. This model has a hierarchical structure which makes it possible to distinguish between target event date θ (the date of interest for the archaeologist) and dates t i of events (artifacts) dated by chronometric methods, typo-chronology or historical documents. One assumes that these artifacts are all contemporaneous, that is relevant to the date of the target event. The dates can be affected by irreducible errors, hence the possible presence of outliers. To take into account these errors, the discrepancy between dates t i and target date θ is modeled by an individual variance σ 2 i , which allows the model to be robust to outliers in the sense that individual variances act as outlier penalization. The posterior distribution of the variance σ 2 i indicates if an observation is an outlier or not. Thanks to this modeling, it is not necessary to discard outliers because the corresponding high values σ 2 i will automatically penalize their contributions to the event date estimation. Moreover, this model does not require additional exogenous or hyperparameters. The only parameter involved in prior shrinkage, s 2 0 , comes uniquely from the data analysis via the individual calibration process. So, the approach is adapted to very different datasets. The good robustness properties of the event date model are paid with less precision in the dates. However, this loss of precision is compensated by better reliability of the chronology. The event model constitutes the basic element in our chronological modeling approach. Dating data are nested within target event dates (with or without stratigraphic constraints between them) which in turn may be nested into phases (with or without succession constraints between them). Succession constraint, maximal duration and/or minimal hiatus can be put on the event dates in the phases.

MCMC computation

.1. Algorithms. The posterior distributions of the parameter of interest θ and of other related parameters t i and σ i can not be obtained explicitly. It is necessary to implement a computational method to approximate the posterior distributions, their quantiles, the Bayes estimates and the highest posterior density (HPD) regions. We adopt a MCMC (Markov Chain Monte Carlo) algorithm known as the Metropolis-within-Gibbs strategy because the full conditionals cannot be simulated by standard random generators. For each parameter, the full conditional distribution is proportional to [START_REF] Vibet | Bayesian Modeling of Archaeological Chronologies[END_REF]. Details on the algorithms used are given in [START_REF] Lanos | Hierarchical Bayesian modeling for combining dates in archaeological context[END_REF].

Here we give a more detailed insight on the way to estimate the date t i which is defined on the set R due to the random effect model chosen in (2). The full conditional distribution of t i (∀i = 1, ..., n) is given by ( 34)

p(t i | q) ∝ 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 2 i (t i ) = s 2 i + σ 2 gi (t i ). Symbol q represents the observations and all other parameters according to equation [START_REF] Vibet | Bayesian Modeling of Archaeological Chronologies[END_REF].

To estimate the posterior density of t i , we can choose between three MCMC algorithms. In each case, the support of the proposal density includes the support R of the target posterior density.

• MH-1: A Metropolis-Hastings algorithm where the proposal is the prior distribution of the parameter. This method is recommended when no calibration is needed, namely for TL/OSL, Gaussian measurements or typochronological references. • MH-2: A Metropolis-Hastings algorithm where the proposal is an adaptive Gaussian random walk. This method is adapted when the density to be approximated is unimodal. The variance of this proposal density is adapted during the process. • MH-3: Metropolis-Hastings algorithm where the proposal mimics the individual calibration density. This method is adapted for multimodal densities, such as calibrated measurements. We are frequently confronted with multimodal target distributions as, for instance, in archaeomagnetic dating (see example 6.1). In this case, algorithms MH1 and MH-2 are not well adapted to ensure a good mixing of the Markov chain. Consequently, an alternative is to choose the MH-3 algorithm with a proposal distribution that mimics the individual calibration density defined in (7). In order to ensure the convergence of MCMC algorithm, the support of the proposal must be R. Thus we can consider a mixture having a gaussian component. This component ensures that the whole support is visited. So we take as proposal :

(35) λC + (1 -λ)N (0, 1 4 (T b -T a ) 2 )
where λ is a number fixed close to 1, and the distribution C approximates the individual calibration density. We can choose for instance the empirical measure calculated on the simulated sample from (7) or a mixture of uniform distribution: Figure 14 and 15 provide graphical tools for the diagnostic of MCMC sampler. We consider the example of Lezoux and only give the result for the parameter of interest θ and the date t dec associated with the declinaison measurement in AM dating. Note that the behavior is the same for all the parameters of the model.

Whithin the Gibbs sampler, the t dec is simulated using MH-3 algorithm, and so the proposal distribution depends on the choice of parameter λ in 35. For three values of λ, Figure 14 and 15 show a stationary behavior with good mixing properties. The autocorrelations even at lag 1 are small enough, to provide a good approximation of the posterior distribution and its characteristics (mean/ variance

Figure 1 .

 1 Figure1. Graphical representation of the densities defined in (12) and (14) with the following parameter values: s j = 30 , σ g = 10, σ δ = 10 2 and θ = t = 1000

Figure 2 .

 2 Figure 2. DAG for the hierarchical model of a group of event dates.

  -T a ) 2 .1 P (α, β)with P = {(α, β) | T a ≤ α ≤ β ≤ T b }.The pairs (α, β) are uniformly distributed on the triangle P .

Figure 3 .

 3 Figure 3. "NSBC" phase model: prior joint density on two (non ordered) event dates θ 1 and θ 2 between two boundaries α and β and after integration against α and β: a concentration effect appears around values θ 1 = θ 2 .

Figure 4 .

 4 Figure 4. Distributions of the span ∆ pour different values of the number r of event dates in a NSBC phase model with: T a ≤ α ≤ θ 1 , . . . , θ r ≤ β ≤ T b . When r is small (r = 2, 3, 4), the densities are very high for ∆ near zero, reflecting a high concentration effect.

Figure 5 .

 5 Figure 5. An other way to look at the concentration effect in the NSBC phase model: distribution of (32) for r = 10 event dates with: T a = 0 ≤ α ≤ θ 1 , . . . , θ 10 ≤ β ≤ 100 = T b .

Figure 6 .

 6 Figure 6. Posterior distribution of parameters in NSBC phase model [left] and in event date model [right].

Figure 7 .Figure 8 .

 78 Figure7. Evolution of the Bayes estimate of θ (the date of the target event (5)) as a function of the proportion q of outlier for three Bayesian models. The number of replication is 500. The data are simulated from with µ = 5 and the sample size is 100. The hyperparameter in (12) are fixed σ δ = 10 and p = 5% .

Figure 9 .

 9 Figure 9. Boxplot of the Bayes estimates of parameters (θ i , t i , σ i ) i=1,...,10 in the Target event model (3). The circles (resp. stars ) indicate the true values of the dates θ i . (resp t i ) The number of replication is equal to 500.

Figure 11 .

 11 Figure 11. Lezoux . [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 θ. The bar above the density represents the shorter 95% posterior probability interval (credibility interval). The vertical lines, delimiting the colored area under the density curve, indicate the endpoints of the 95% highest posterior density (HPD) region.

Figure 12 .

 12 Figure 12. Lezoux. Posterior densities obtained for standard deviations σ i .

Figure 13 .

 13 Figure 13. Malpaís Prieto, without informative constraint [left] and with constraint of maximal duration τ 0 = 60 [left]. Each picture provides the posterior densities for event dates θ j [gray background] and the posterior densities for beginning and end of the phase [dark gray background]. The bar above the density represents the shorter 95% posterior probability interval (credibility interval).

  ) i=1,...,M are the ordered values of the simulated sample. An alternative is to choose the distribution with density1 C p i=1 p i (τ i )I [τi,τi+1] ,whereC = p i=1 f (τ i )(τ i+1 -τ i )where (τ i ) i is a deterministic grid of T and p i is the density of individual calibration density (7). .2. Assessing Convergence. Using the packages ArchaeoPhases and coda implemented in R (see Philippe and Vibet (2017b); Plummer et al. (2006) and R Core Team (2017)), we check the convergence of MCMC samples simulated for different values of the parameters λ in (35).

Table 1 .

 1 95% HPD region for estimating the phase (begin/end/duration)

	PHASE
	Begin	[1209; 1355]
	End	[1372; 1450]
	Duration	[48; 211]

Table 2 .

 2 95% HPD regions for estimating the phase (begin/end/duration)

	, right) :

Acknowledgments

The authors expressed their gratitude to Grégory Pereira (CNRS, UMR ArchAm, University of Paris 1 Panthéon-Sorbonne) for the exchange of information about the example of "Malpaís Prieto site in Mexico". We thank Marie Anne Vibet for her very helpful discussions during the preparation of this paper.

/ quantiles etc) . Morever it is important to note that the results are not sensitive to the value of λ.

The Gelman diagnostic (evaluated from 5 parallel chains) is equal two 1, and so it confirms the convergence of MCMC samplers.