An Archaeological Problem

As a first example, consider an MCMC analysis of the model and data described in Naylor and Smith (1988). This concerns inference about five dates forming boundaries between differing, distinct, and abutting periods of pottery production. The data consist of radiocarbon date determinations associated with pottery fragments identified as belonging to one of four periods of pottery production.

Their analysis obtains a posterior density


\begin{displaymath}
p({\mbox{\boldmath$\alpha$}} \vert {\mbox{\boldmath$x$}})=\...
...$}})}p({\mbox{\boldmath$\alpha$}})}{p({\mbox{\boldmath$x$}})}
\end{displaymath} (10)

where

${\mbox{\boldmath$\alpha$}}$ is a vector of dates before present2 (BP) $ \alpha_1 >
\alpha_2 > \alpha_3 > \alpha_4 > \alpha_5$, having prior


$\displaystyle p({\mbox{\boldmath$\alpha$}})$ $\textstyle =$ $\displaystyle 1~~~\mbox{ if }
\alpha_1 > \alpha_2 > \alpha_3 > \alpha_4 > \alpha_5 > 0$ (11)
  $\textstyle =$ $\displaystyle 0~~~\mbox{ otherwise }$ (12)

$p({\mbox{\boldmath$x$}})$ is the normalizing constant.

$L({\mbox{\boldmath$\alpha$}} ; {\mbox{\boldmath$x$}})$ is the Log-Likelihood


\begin{displaymath}
\sum_{i=1}^n \log p(x_i\vert s_i,j_i, {\mbox{\boldmath$\alpha$}})
\end{displaymath} (13)

where


$\displaystyle p(x_i\vert s_i,j_i, {\mbox{\boldmath$\alpha$}})$ $\textstyle =$ $\displaystyle \int_{\alpha_{j+1}}^{\alpha_j}
p(x \vert \mu (\theta ), s)p(\theta\vert \alpha_j, \alpha_{j+1})d\theta$ (14)
  $\textstyle =$ $\displaystyle (\alpha_j - \alpha_{j+1})^{-1} \sum_{i=l_i(j)}^{l_2(j)+1}I_i$ (15)
$\displaystyle \mbox{with}$     (16)
$\displaystyle I_i$ $\textstyle =$ $\displaystyle \int_{d_{i-1}(j)}^{d_i(j)} p(x\vert a_i + b_i\theta,s)d\theta$ (17)
  $\textstyle =$ $\displaystyle b_i^{-1}\left\{\Phi\left[\frac{a_i+b_id_i(j)-x}{s}\right]
-\Phi\left[\frac{a_i+b_id_{i-1}(j)-x}{s}\right]\right\}$ (18)

$\Phi(\cdot)$ being the standard Normal distribution function and $d_i, a_i, b_i, l_i(j) \mbox{ and } \mu(\theta)$ are values expressing a piecewise linear radiocarbon calibration curve (Clark, 1979).


$\displaystyle p(\theta\vert j, {\mbox{\boldmath$\alpha$}})$ $\textstyle =$ $\displaystyle p(\theta \vert \alpha_j, \alpha_{j+1} ) ~~~
\alpha_j > \theta \geq \alpha_{j+1}$ (19)
  $\textstyle =$ $\displaystyle 0~~~\mbox{ otherwise }$ (20)

where here $\theta$ is the actual date of manufacture.

The analysis presented here, using MCMC, incorporates into the likelihood a calibration curve, defined by Clark (1979) in piecewise linear segments, relating radiocarbon date to calendar date. This curve includes inversion (i.e. it is not monotonically increasing), which renders its use difficult. Figure 1 shows univariate posterior densities estimated using BKDE for four of the five boundary dates.

Taking ${\mbox{\boldmath$\theta$}} = \log(h)$ (ensuring that the estimate of $h$ at any point is positive) to have prior $N(0,1)$, represents a loosely held belief ($\sigma = 1$) that the density is smooth ($\mu
= 0$), corresponding to $h=1$. That is $h$ in the range $0.135$ to $7.389$ being considered probable a priori.

These densities are seen to be considerably less smooth than the corresponding results of Naylor and Smith (1988) who used an older, but smoother, calibration curve. There is some evidence that the likelihood is not smooth and other authors have found that there are considerable problems for maximum likelihood analyses (Clark, 1979; Cunliffe, 1984).

The analysis returns as many points as required. In this example $10,000$ points were generated and then reduced to 500 by extracting every 20th result 3. Taking $n_0=1$, a uni-dimensional likelihood is constructed with one previous sample.

Figure:Archaeological problem, Bayesian estimation.
\begin{figure}
\centering
\psfig{figure=/home/danny/thesis/pics/arc.ps,width=5.25in,angle=270}
\end{figure}

Figure:Archaeological problem, different priors for ${\mbox{\boldmath $\theta$}} = \log(h)$.
\begin{figure}
\centering
\psfig{figure=/home/danny/thesis/pics/two.ps,width=5.25in,angle=270}
\end{figure}

Figure:1998 International $^{14}C$ atmospheric data set (24,000 to 0 BP).
\begin{figure}
\centering
\psfig{figure=/home/danny/thesis/pics/c14fig.ps,width=5.25in,angle=270}
\end{figure}

In Figure 2 the effect of changing the prior is observable:

  1. Prior $N(0,1)$ loosely held belief in a smooth density, ( $0.135~\le~h~\le~7.389$).
  2. Prior $N(0.1,1)$ loosely held belief in a slightly smoother density ( $0.149~\le~h~\le~8.166$).
  3. Prior $N(0,0.1)$ stronger belief in a smooth density ( $0.729~\le~h~\le~1.372$).
  4. Prior $N(0.1,0.1)$ stronger belief in a slightly smoother density ( $0.806~\le~h~\le~1.516$).

On initially examining Figure 2, it is apparent that 1 and 2 show the same curve, the weak prior is overridden by the data. In 3 and 4 the curves are similar and smooth, the information in the data being insufficient to override the strong prior.

Having a degree of belief in the confidence placed in an ``expert opinion'' is central to a subjective Bayesian approach to a problem. Priors arising from either previous knowledge of a process or from some reliable analysis of a situation may be assumed to have some degree of validity. This is expressed here in 3 and 4 by the adoption of a small value for $\sigma$ leading to the retention of the smooth curve in the density estimate. The larger prior value for $\sigma$ in 1 and 2 expresses the opposite, i.e. an assumption that there is no prior belief in a smooth curve and that is reflected in the density estimation obtained. Equally a strong belief in a jagged curve could be imposed, which would seem to be more appropriate here. In all cases, more data reduces the influence of the prior.

Finally, note that the jagged curve in 1 and 2 might be taken to indicate that the curve is under smoothed, however, the $^{14}C$ calibration curve used here is in nowise linear as can be seen in Figure 3, and this might be reflected in the posterior.

danny 2009-07-23