[Back] [Forward] [Contents] [Journal Homepage]

4 Bayesian wiggle-matching (technical)

This section is intended for statisticians.

The theory of Bayesian wiggle-matching has been developed in full in Christen and Litton (1995), and in general, the theory of Bayesian Calibration has been developed in Buck et al. (1991); (1992); Christen (1994a); Christen and Buck (1998); Buck and Christen (1998), Gomez-Portugal Aguilar et al. (2002); Nicholls and Jones (2001), among others. Here the theory of wiggle matching is briefly presented, plus a further consideration to deal with outliers.

Let y = (y1, y2, ym) be a series of radiocarbon determinations with their corresponding standard errors σ1, σ2, ... ,σm, and their associated calendar ages θ = (θ1, θ2, ... , θm) (that is, the calendar ages at which each sample radiocarbon dated ceased metabolising). The output of a radiocarbon analysis, as provided by the laboratory, is an estimated date in the radiocarbon scale (yj) and a standard error (σj). The standard error reported is calculated using both empirical and theoretical considerations and the usual assumption is to consider it as known (see below). Here, we assume that θi - θj it is known exactly; most likely the radiocarbon dated samples were taken from chunks of tree-rings in a log or timber. The model we use is the following

LaTex math formula:
$$
y_i \mid \theta_j \sim N\left(
\mu(\theta_j), \sqrt{\sigma^2_j + \sigma^2(\theta_j)}\right),
$$

where μ (θj) is the piece-wise linear calibration curve (we are using the INTCAL98 calibration curve, see Stuiver et al. 1998) and σ2j) is a variance term arising from the errors observed in the calibration curve (see Christen and Litton 1995 and Christen 1994a). We also assume that ƒ(y | θ) = π mj=i ƒ(yi | θj) (conditional independence). That is, σj is assumed known, which is the usual assumption in radiocarbon dating, and the determinations are considered as conditional independent.

The piece-wise linear calibration curve is defined as

LaTex math formula:
$$
\mu(\theta) =
\left(\frac{\theta - t_{k-1}}{t_k - t_{k-1}}\right)x_k +
\left(\frac{t_k - \theta}{t_k - t_{k-1}}\right)x_{k-1}
$$

and the variance term is defined as

LaTex math formula:
$$
\sigma^2(\theta) =
\left(\frac{\theta - t_{k-1}}{t_k - t_{k-1}}\right)^2\sigma^2_k +
\left(\frac{t_k - \theta}{t_k - t_{k-1}}\right)^2\sigma^2_{k-1} +
\frac{\lambda^2 (\theta - t_{k-1})(t_k - \theta)}{t_k - t_{k-1}},
$$

for tk > θ >= tk-1; where tk ± σk is the pooled radiocarbon determination for the calibration curve at knot (calendar year) xk. λ is taken to be 19 (see Christen and Litton 1995 and Christen 1994a), although Gomez-Portugal Aguilar et al. (2000) argue that a lower value could be used.

The definition of σ 2(θ) is based on a Brownian bridge model (see Billingsley 1999, 93) and represents a simple, although realistic, way to include the errors in the calibration curve (see Christen and Litton 1995; Christen 1994a; Gomez-Portugal Aguilar et al. 2002 and Nicholls and Christen 2000 for a more in-depth discussion). The option where no errors are included in the curve is covered simply by fixing σ 2(θ) = 0.

In many applications, and for practical reasons, the inclusion of the errors in the calibration curve is neglected. This is in fact good practice, as the variance in the posterior is comparatively higher than the variances in the curve. This is normally not the case in wiggle-matching because of the increased precision arising from the simultaneous calibration of several radiocarbon dates, even if the individual standard errors are big.

Let us assume that θ1 <= θ2 <= ... θm and that θ0 is the (unknown) calendar date for the event to be dated (this event may or may not coincide with any of the θj). We assume that θj = θ0 + αj where, as mentioned above, αj is known exactly (αj could be negative). The only unknown parameter then is θ0 and its posterior density is given by

LaTex math formula:
$$
f( \theta_0 \mid \vy ) \propto
f(\theta_0) \prod_{j=1}^m \frac{1}{w_j(\theta_0 + \alpha_j)}
e^{-\frac{\left\{y_j - \mu (\theta_0 + \alpha_j)\right\}^2}
	  {2w_j^2(\theta_0 + \alpha_j)}}
$$

where wj2 (θ) = σj2 + σ 2(θ). The prior f(θ0) is an uniform distribution provided by the user (upper and lower bounds for θ0). This posterior is normalised numerically and provided as a histogram, both with (see Figure 2) and without (see Figure 1) considering the errors in the calibration curve. (Unfortunately Bwigg only allows for uniform priors, although mexcal could be used directly to consider more sophisticated priors. This feature is not yet included in Bwigg.)

Outliers

The analysis of 'shift outliers' is based on Christen (1994b). A shift in the radiocarbon scale δj is considered for each determination, having a priori distribution

LaTex math formula:
$$
\delta_j \sim N( 0, \beta \sigma_j^2)
$$

We fix β = 2. The latent variable Φj = 0,1 is also introduced. If yj is an outlier then Φj = 1, and thus needs a shift δj in the radiocarbon scale to be corrected. The following model [2] is assumed:

LaTex math formula:
$$
y_i \mid \theta_j,\delta_j, \phi_j  \sim
N( \mu(\theta_j) + \delta_j \phi_j, w_j(\theta_j))
$$

Interest then focuses on approximating P(Φj = 1 | y), the posterior probability of determination j being an outlier. MCMC is used to approximate P(Φj = 1 | y). In turn, a simulation is run from the full conditionals of each of the Φj's (Bernoulli distributions) and δj's (normal distributions, see Christen 1994b). To simulate from θ0, a proposal is simulated θ0' form f(θ0) and accept it with probability

LaTex math formula:
$$ \min \left\{ 1, \frac
{
\prod_{j=1}^m \frac{1}{w_j(\theta_0' + \alpha_j)}
e^{-\frac{\left\{y_j - \mu (\theta_0' + \alpha_j)\right\}^2 - \delta_j \phi_j}
	  {2w_j^2(\theta_0' + \alpha_j)}}
}
{
\prod_{j=1}^m \frac{1}{w_j(\theta_0 + \alpha_j)}
e^{-\frac{\left\{y_j - \mu (\theta_0 + \alpha_j)\right\}^2 - \delta_j \phi_j}
	  {2w_j^2(\theta_0 + \alpha_j)}}
}
\right\}
$$

The resulting MCMC is very well behaved in that the support of f(θ0) is not too wide in comparison to the main bulk of the posterior for θ0. 10,000 samples are taken, every 10 passes, with a burn-in of 1,000 passes. This is usually enough for approximating P(Φj = 1 | y).


Footnote 2: Note that the errors in the calibration curve are taken into consideration via wjj)


[Back] [Forward] [Contents] [Journal Homepage]

© Internet Archaeology URL: http://intarch.ac.uk/journal/issue13/2/4tech.html
Last updated: Thu Jan 23 2003