# Inference for population dynamics in the Neolithic period\thanksrefT1

###### Abstract

We consider parameter estimation for the spread of the Neolithic incipient farming across Europe using radiocarbon dates. We model the arrival time of farming at radiocarbon-dated, early Neolithic sites by a numerical solution to an advancing wavefront. We allow for (technical) uncertainty in the radiocarbon data, lack-of-fit of the deterministic model and use a Gaussian process to smooth spatial deviations from the model. Inference for the parameters in the wavefront model is complicated by the computational cost required to produce a single numerical solution. We therefore employ Gaussian process emulators for the arrival time of the advancing wavefront at each radiocarbon-dated site. We validate our model using predictive simulations.

10.1214/12-AOAS579 \volume6 \issue4 2012 \firstpage1352 \lastpage1376

Population dynamics in the Neolithic

T1Supported in part by Leverhulme Grant F/00125/AD.

A]\fnmsAndrew W. \snmBaggaleylabel=e1]a.w.baggaley@ncl.ac.uk, A]\fnmsRichard J. \snmBoys\correflabel=e2]richard.boys@ncl.ac.uk, A]\fnmsAndrew \snmGolightlylabel=e3]andrew.golightly@ncl.ac.uk, A]\fnmsGraeme R. \snmSarsonlabel=e4]graeme.sarson@ncl.ac.uk and A]\fnmsAnvar \snmShukurovlabel=e5]anvar.shukurov@ncl.ac.uk

Bayesian inference \kwdpopulation dynamics \kwdthe Neolithic \kwdGaussian process emulation \kwdMarkov chain Monte Carlo

## 1 Introduction

The Neolithic (10,000–4500 BC), the last period of the Stone Age, was the dawn of a new era in the development of the human race. Its defining innovation, the transition from foraging to food production based on domesticated cereals and animals, was one of the most important steps made by humanity in moving toward modern societies. Among dramatic changes resulting from the advent of farming were sedentary living, rapid population growth, gradual emergence of urban conglomerates, division of labor and the development of complex social structures. Another important innovation of the period was pottery making. The mechanism of the spread of the Neolithic remains an important and fascinating question. The most striking large-scale feature of this process is its regular character, whereby farming technologies spread across Europe at a well-defined average speed of about kmyear. Impressive and convincing evidence for this emerged as soon as radiocarbon (C) dates for early Neolithic sites became available [Ammerman and Cavalli-Sforza (1971)], and later studies have confirmed the early results [Gkiasta et al. (2003)]. Figure 1 displays the C dates used in our analysis, and they clearly show the gradual, regular spread of the Neolithic over a time span of over 4000 years from its origin in the Near East around 9000 years ago, to the north-west of Europe.

It is worth emphasising that “the Neolithic” is not a single phenomenon and, although it is often characterized by archaeologists as a “package” of related traits, the individual elements of the package need not have been transmitted simultaneously. This might suggest the need to model each distinct element separately. However, several elements do appear to have traveled together during the spread of farming into most of Europe [Burger and Thomas (2011)] and the local transition to the full package may have been more rapid than has often been assumed [Rowley-Conwy (2011)]. Therefore, we choose to model the Neolithic spread as a single phenomenon.

Despite the overall regular character of the expansion, there are notable regional variations in the spread of the Neolithic [Gkiasta et al. (2003), Bocquet-Appel et al. (2009)]. Careful analysis of the radiometric and archaeological evidence has identified various major perturbations: First, slowing down at topographic obstacles such as mountain ranges. Second, latitudinal retardation above 54N latitude [Ammerman and Cavalli-Sforza (1971)], presumably due to competition with the pre-existing Mesolithic population, whose density was especially high in the north [Isern and Fort (2010)]; the negative influence of the soil type and harsher climate in the north also cannot be excluded. Third, significant acceleration along the Danube–Rhine corridor and the west Mediterranean [Ammerman and Cavalli-Sforza (1971)]; see also Davison et al. (2006). These variations in the rate of spread complicate quantitative interpretations of the C data.

Moreover, the data themselves are far from being complete and accurate. There are inherent errors from the radiocarbon measurement in the laboratory [Scott, Cook and Naysmith (2007)] and uncertainties in the conversion of the C into the calendar age, known as the calibration of the C dates [Buck et al. (1991), Blackwell and Buck (2008)]. Further errors result from archaeological factors, such as uncertain attribution of the dated object to the early farming activity, disturbed stratigraphy of archaeological sites, etc. [Aitken (1990)]. Perhaps even more importantly, objects related to the first appearance of farming at a given location are not necessarily related to the first arrival of farming to its larger region. The occupation of many Neolithic sites could be a secondary process behind the wave of advance. The only obvious method to identify such sites is by continuity with its neighbors.

Given the obviously random nature of the demographic processes underlying the Neolithic expansion, and the abundance of both random and systematic errors in the quantitative (mainly radiometric) evidence available, it is clear that serious statistical analysis is required to quantify the spread of the Neolithic. All previous work in this direction has determined the parameter selection for various models and assessed goodness of fit by using parameter scans. To the best of our knowledge, the present paper is the first to use formal statistical inference methods to determine appropriate ranges for model parameters and assess the fit of the model.

We use a carefully selected set of C dates for the earliest Neolithic sites across Europe to fit a model of the Neolithic expansion. The model relies on the fundamental fact, established empirically as briefly explained above, that the incipient farming spread at a nearly constant average speed. We also include regional variations in the speed associated with the latitudinal gradient, topography and major waterways. Our statistical model adopts an error structure which accounts for both measurement error in the radiocarbon dating process and a spatial error process which smooths departures from the wavefront model. We describe how to obtain the Bayesian posterior distribution for the unknown model parameters using simulations from the deterministic model. Simulations of the expanding Neolithic front are computationally expensive, even at a prescribed deterministic velocity, precluding their use for Bayesian inference. Therefore, we construct a Gaussian process emulator [Santner, Williams and Notz (2003)] for the model, that is, a stochastic approximation to the arrival times obtained from this model. Such methods are widely used in the computer models literature; see, for example, Kennedy and O’Hagan (2001) and references therein. Rather than build a complex space–time emulator to describe the wavefront model, we pragmatically build individual emulators for each site and rely on the spatial Gaussian error process to smooth across geographical space. Finally, we fit this emulator model to our data and describe the inferences we draw.

The remainder of this paper is organized as follows. We describe the C data used in our analysis in Section 2. The implementation of the wavefront model and the statistical model which links the wavefront model to the data are presented in Section 3. The Gaussian process emulator model and its use in constructing an inference scheme are the subject of Sections 4 and 5. Finally, the results of fitting our model to the radiocarbon data are reported in Section 6, and we conclude with a discussion in Section 7.

## 2 Selection and analysis of radiocarbon dates

It is clear from the discussion in Section 1 that C data need to be carefully selected for model parameter inference. Databases of radiocarbon dates contain several thousand entries for Neolithic sites across Europe [e.g., Gkiasta et al. (2003)], but a relative minority of them are suitable for the studies of the initial spread of the Neolithic. In this paper, we use carefully selected C dates from 302 earliest Neolithic sites in Southern and Western Europe from the compilation of Davison et al. (2009), where a detailed discussion of the selection criteria can be found. The objects dated include grains and seeds, pottery, bone, shells and animal horns, wood, charcoal and peat, and soil. We note that, in general, there may be issues in using such a heterogeneous data set, such as those relating to “chronometric hygiene” and the problems that may arise with certain materials or treatments; see, for example, Pettitt et al. (2003) and the discussion in Steele (2010). However, we are not aware of any issues with the provenance of the samples we use in our analysis; see Davison et al. (2009) for further details.

The putative origin of the expansion is in an extended region in the Near East, and, following the results of Ammerman and Cavalli-Sforza (1971), later confirmed by many authors, it is reasonable to place it near Jericho. Similarly to Davison et al. (2009), the starting date and position are adopted as 6572 BC and N, E), the Tell Kashkashok site near Jericho.

Many sites in the selection have several dates (often 3–10, and occasionally 30–50) which can be treated as the arrival time contaminated by errors. Each date is published together with the accuracy of radiocarbon measurement in the laboratory, but we recall that other sources of random and systematic errors are present as well. Suppose that, at site , we have calibrated C dates (, ) and their standard deviations . These data are calculated from the central 99.7% interval of the calibrated probability distributions for each object, obtained using calibration curve IntCal04: Reimer et al. (2004). Specifically, each is taken to be the midpoint of its calibrated interval and to be one-sixth of the length of the interval. This data set can be found in the supplementary material [Baggaley et al. (2013)] for this paper. We can calculate summaries for the data at each site, namely, an arrival time and its accuracy , as defined in this context by Scott, Cook and Naysmith (2007), by appropriately weighting the data, as

(1) |

Clearly, there is some loss of information in not using the full calibrated probability distribution for each object. However, to include this within our analysis by using, for example, a normal mixture model representation for each distribution, would make our analysis considerably more complicated. In any case, the model we do use in our analysis is one for the weighted means , whose calibrated distributions will be roughly normally distributed due to the central limit theorem. The left-hand panel of Figure 1 shows the sites on the map, color-coded according to the arrival time .

## 3 Justification and implementation of the wavefront model

The wave of advance model^{2}^{2}2To avoid confusion, we note that
the word “wave” does not refer here to any oscillatory behavior.
has a deep mathematical basis. Such fronts are the salient feature of
a wide class of phenomena that involve an exponential growth of some
quantity (in this case, the population density) and its spread via
random walk or related diffusion; see Fort and Pujol (2008) and references
therein. This behavior is captured by the classical equation of
Fisher (1937) and Kolmogorov, Petrovskii and
Piskunov (1937), known as the FKPP equation,
of the form

(2) |

where is the population density, a function of position and time , is the growth rate of the population (the difference of the birth and death rates per unit area), is the carrying capacity of the landscape (the maximum sustainable population density), and is the diffusivity, a measure of human mobility. For constant , this equation has obvious solutions and . Solutions of the initial value problem for this equation in a homogeneous domain, with a localized initial condition, have the form of a wave of advance, in one dimension, where the solution changes from ahead of the wavefront to behind it. The width of the transition region (an internal boundary layer) is of order and, in one dimension, the front propagates at a constant speed

(3) |

independently of . In two dimensions, the constant-speed propagation provides a good approximation as soon as the radius of curvature of the front becomes much larger than . For parameter values characteristic of the Neolithic expansion—for example, year and km s (see the discussions in Sections 3.2 and 3.4), providing km s—the front width is km, and the model of a wavefront propagating at a constant speed is safely applicable at the continental scales of 100–1000 km. Even though we aim at an inference for the speed of the Neolithic expansion, we find it convenient to present our results with reference to the FKPP model. In particular, is parameterized in terms of and since these quantities admit a direct interpretation in terms of human behavior. Nevertheless, we emphasize that the inference that follows does not depend in any fundamental way on the FKPP equation, or on the assumptions of the demographic processes underlying it; the results we present may simply be considered as providing statistical constraints on the speed of a wavefront , independently of any interpretation of the mechanisms behind this wave.

### 3.1 The wavefront representation

A discrete version of the wavefront model, suitable for numerical simulations, is constructed as follows. We select the point source of the population and then, at a small radius from this point, we arrange a small number of test points (“particles”) on a circle around it, with the position of particle denoted by , where and are the geographical latitude and longitude, respectively. This allows us to compute the local tangent and normal vectors of the front, in particular, the local unit outward normal (in the direction of the advancing front). At each time step we move each particle along with the local velocity given by (6). The wavefront thus expands outward in time; in a homogeneous, isotropic domain, the expansion would take the form of concentric circles, but in the heterogeneous and anisotropic domain used here (described in detail below), the expansion is irregular and tracks the front obtained from the corresponding FKPP model. At the continental scale, the Earth’s curvature is noticeable and we work on a spherical surface of the appropriate radius.

The separation of the particles increases as the front expands, and an additional particle is introduced between any two particles if their separation exceeds a specified value . Conversely, if the front contracts, particles that are closer than are removed. We take to be equal to the grid spacing (introduced below) of 4 arc-minutes, that is, between 3 km at 60N and 7 km at 30N. It is important to monitor the particle ordering along the front after adding or removing particles in order to maintain a continuous wavefront.

When the front encounters an obstacle (e.g., a mountain range), its flanks move round it and then have to merge behind the obstacle. The merger is implemented by applying algorithms initially used to model the evolution of magnetic flux tubes in astrophysical simulations, and so described in more detail in Baggaley et al. (2009). Briefly, we check, at each time-step, if any two particles (which are not neighbors) are closer than ; if so, the particle ordering is switched to “short-circuit” any (almost closed) loops that have arisen. The front then propagates further, leaving behind closed loops which are then removed from the simulation since they do not affect the front arrival time.

### 3.2 Regional variations

The front speed (3) depends on the product of and , so that these parameters cannot be estimated simultaneously from studies of the front propagation. The growth rate is rather well constrained by biological and anthropological factors: Birdsell (1957) reports values in the range 0.029–0.035 year and Steele, Adams and Sluckin (1998) suggest the range 0.003–0.03 year. Here, for comparability with our own recent studies [e.g., Davison et al. (2006)], we take year, corresponding to a generation time of 30 years. This choice does not limit the applicability of our model. If, for example, a slightly different value of is preferred, then the reported values of need only to be scaled so that the product remains unchanged.

To make the problem tractable, we adopt a specific spatial variation of the diffusivity and seek inference for its magnitude, just a single scalar parameter. Such assumptions unavoidably involve significant arbitrariness; however, a well-informed assumption is justifiable if it results in a model that fits the observations to a satisfactory degree. Thus, we represent the diffusivity in the form

(4) |

where is the constant-dimensional magnitude and is a dimensionless function of position whose magnitude is of order unity and whose form is assumed to be known. We prescribe to depend solely on the topography (altitude above the sea level) and latitude as described below.

Early farming does not appear to have been practical at high altitudes in the temperate zones considered here [e.g., sites in the Alpine foreland are not found on land above 1 km: Whittle (1996)]; to include this effect in our model, we smoothly reduce to zero at around this height. Also, in order to allow for coastal sea travel, exponentially decreases with the distance to the nearest land. Finally, we include a linear decrease of toward the northern latitudes. Altogether, we adopt the following form for the dimensionless diffusivity:

(5) |

where is the altitude, is the latitude, and we recall that corresponds to the seas. Both and are complicated functions of position reflecting the topography and coastlines of Europe. Of course, the choice of these particular forms are only a crude and exploratory attempt to capture the dependence on latitude and altitude. However, we have found our results not to be too sensitive to minor changes in these functions.

The altitude data have been taken from the ETOPO1 1-minute Global Relief database [Geophysical Data System (2011)]. As a reasonable compromise between computational efficiency and spatial accuracy, we use a grid with a spatial resolution of 4 arc-minutes between W and E in longitude and N and N in latitude. The boundaries of the domain explored and the spatial variation of are shown in Figure 1.

Significant regional anomalies revealed by archaeological and radiometric evidence are the early Neolithic Linearbandkeramik (LBK, Linear Pottery) Culture that propagated in 5500–4900 BC along the Danube and Rhine rivers at a speed of at least 4 kmyear [Dolukhanov et al. (2005)] and the Impressed Ware culture (5000–3500 BC) that spread along the west Mediterranean coastline at an even higher speed. The latter could be as high as 10 kmyear [Zilhão (2001)]. Davison et al. (2006) included directed spread (advection) along the major river paths and coastlines into their mathematical model to achieve a significant improvement in the fit to C dates.

We include advection along the Danube–Rhine river system and the coastlines using world map data from Pape (2004), in the form of tangent vectors of the coastlines and river paths, and , respectively, given at irregularly spaced positions. To remap the data to the regular grid introduced above, we use the original data weighted with , where is the distance between the grid point and the data point, and the chosen scale of 15 km contains 2–5 grid separations at latitudes 30–60N.

With the unit tangent vectors for the coastlines and river paths thus defined, the magnitudes of the respective local front velocities, and , are subject to inference. In other words, the local front velocity is given by

(6) |

where . Since, unlike the geographical data, the points defining the front are not necessarily at the nodes of the regular grid, we used bilinear interpolation from the four closest grid points to calculate (and thus ), and at the front positions. We note that the accuracy of the wavefront model was tested by comparing it with numerical solutions of (2) for a wide range of parameter values; the agreement was invariably quite satisfactory with a typical deviation of 40 years.

### 3.3 Statistical model

We model the radiocarbon dates for each object at each site as being a noisy version of those dates predicted by the wavefront model, explicitly allowing for three types of error introduced in Section 1: the date of object at a site (location ), for , , is modeled as

where is the deterministic wavefront arrival time, a function of , and and are independent error terms following a standard normal distribution. Here allows for spatially inconsistent data, is the standard deviation of the C measurement in the laboratory and calibration, and is a global error term. Between them, the and terms allow for mismatches between our model and the data: allows for a certain level of variability between the arrival times at nearby sites, while allows for an additional global variability (with no spatial correlation). As well as the mismatch arising when a simple large-scale model is applied to a process with inherent variability on smaller scales, these terms will also allow our model to be fairly robust to any problems in our radiocarbon data such as the misinterpretation of sites of secondary Neolithic occupation as being indicative of the wave of first arrival. We note that some authors [e.g., Buck et al. (1991)] choose to model the start-date of a “first arrival phase” using an additional parameter for each site, with the explicit expectation that all individual samples will fall after this date. However, our wavefront model attempts to capture the Neolithic transition over much larger spatial and temporal scales and so we choose not to incorporate such detailed temporal effects.

We model the spatially smooth error process using a zero mean Gaussian process (GP) prior [Rasmussen and Williams (2006)] with covariance function , that is,

We impose smoothness by taking a Gaussian kernel for the covariance function so that the covariance between spatial errors at locations and is

The parameters of this function control the overall level of variability and smoothness of the process, with larger values of the length scale giving smoother realizations and thereby down-weighting sites with different dates relative to nearby sites.

In the following analysis it will be more straightforward to work with an equivalent model for the summary statistics in (1), namely,

(7) |

where the and are independent error terms following a standard normal distribution. Thus, the inferential task is to make plausible statements about the wavefront model parameters , the global error term and the Gaussian process hyperparameters and .

### 3.4 Bayesian inference

We adopt a Bayesian approach to inference and express our prior uncertainty for the unknown quantities through a density . The speed of spread reported by Ammerman and Cavalli-Sforza (1971), together with the growth rate discussed in Section 3.2, suggests kmyear. The spread reported for the LBK culture [e.g., Ammerman and Cavalli-Sforza (1973), Gkiasta et al. (2003)] suggests kmyear. The spread reported for the Impressed Ware culture [Zilhão (2001)] would suggest kmyear. However, this study is restricted to the Western Mediterranean and we believe this parameter would take a smaller value over the European continent such as kmyear. We do not expect the wavefront model to fit the data well in all regions, so we expect a relatively large value of (compared with the accuracy of individual dates), years. The magnitude of , which models the variability in dates between spatially nearby sites, is difficult to estimate, but the maximum accuracy of laboratory radiocarbon measurements (which will probably be the smallest source of variability within our model) suggests a minimum significant value of years. The mean minimum separation of radiocarbon-dated sites suggests . We use the values above to guide our choice of mode for the prior distribution. In view of the uncertainties in these values, we make the prior rather diffuse, with independent components. Specifically, we use the log-normal and inverse-gamma distributions, with

Let and denote the observed arrival times and those from the wavefront model (evaluated at ). Then, using (7), the data model is an -dimensional normal distribution, with

where , has th element , and are the locations of the radiocarbon-dated sites. Combining both sources of information (the data/model and the prior) gives the joint posterior density as

In practice, the posterior density is analytically intractable and we therefore turn to a sampling-based approach to make inferences. Markov chain Monte Carlo (MCMC) methods can readily be applied to this problem, but such schemes will typically need many thousands or millions of evaluations of , each requiring a full simulation of the wavefront model expanding across the whole of Europe. This takes around 10 seconds on a quad-core CPU with a clock speed of 2.67 GHz and the code parallelized using OpenMP. This computational cost prohibits the use of an MCMC scheme for sampling from the posterior distribution of the unknown quantities. To proceed, we therefore seek a faster approximation of the first arrival times from the wavefront model.

## 4 An approximation to the wavefront model

The first arrival times from the wavefront model can be approximated by using both deterministic and stochastic methods such as cubic splines and Gaussian processes, respectively. These approximations (known as emulators) are determined by first setting up an experimental design consisting of a number of locations and parameter values, and then computing the (deterministic) first arrival times from the wavefront model using the locations and parameter values in the design. Finally, the proposed emulator is fitted to this model output. In line with many other authors, for example, Kennedy and O’Hagan (2001) and Henderson et al. (2009), we favor using Gaussian processes to emulate the first arrival times, as they not only fit this model output exactly and interpolate smoothly between design points but also quantify levels of uncertainty around interpolated values. Rather than build a complex emulator whose inputs are both site location and wavefront model parameters, we pragmatically build individual (parameter only) emulators for each site and rely on a spatial Gaussian error process to smooth the resulting spatial inhomogeneities. This strategy is viable because the inference scheme described later in Section 5 only requires emulation at the 302 sites in the radiocarbon data set. One possible drawback is that it does produce site-specific emulators that are not as spatially smooth as the more complex emulator. However, there are significant computational benefits of building separate emulators for each of the 302 sites, as they work in a smaller input space and can be computed in parallel.

### 4.1 A Gaussian process emulator

The output from a single simulation of the wavefront model with input parameter values gives the front arrival time for all sites. Consider the arrival time at a single site . We model the emulator for the arrival time at this site using a Gaussian process with mean function and covariance function , that is,

(8) |

A suitable form for the mean function can be determined by noting that the great-circle distance from the source of the wavefront is approximated by , where is the front speed given by (6) and is the arrival time. Hence, and therefore we seek

where are coefficients to be determined, which account for the relative importance of diffusivity, coastal and river speeds. We determine them using ordinary least squares fits. Note that the mean function depends on to be consistent with (3). We use a stationary Gaussian covariance function

where is the length scale at site along the -axis, as is widely used in similar applications since it leads to realizations that vary smoothly over the input space [Santner, Williams and Notz (2003)]. We describe how to make inferences about the parameters of these site-specific covariance functions in Section .2.

Suppose that simulations of the (computationally expensive) wavefront model are available to us, each providing the arrival time at each radiocarbon-dated site. Let denote the -vector of arrival times resulting from the wavefront model with input values , where . From (8) we have

where is the mean vector with th element and is the variance matrix with th element .

We can use the simulations to model the front arrival times corresponding to other values of the input parameters. Indeed, suppose we need to simulate front arrival times at a collection of new design points . This is straightforward as, with the local emulators, the arrival times are mutually independent and the distribution of the arrival time at site conditional on the training data can be determined using standard properties of the multivariate normal distribution as

(9) |

where

(10) | |||||

(11) |

Note that, to simplify the notation, we have dropped the dependence in these expressions on the hyperparameters , where .

## 5 Inference

Before we can fit the Gaussian process emulators, we need to obtain the training data from the wavefront model. These data are obtained by running the emulators at a particular experimental design; see Section .1 for details. Given the C arrival times and the training data, it is possible in theory to fit the emulator and statistical model (7) jointly using an MCMC scheme. However, this is likely to be extremely computationally expensive, and we therefore fit the emulator and the statistical model separately. This approach is advocated by Bayarri et al. (2007) and adopted by Henderson et al. (2009), among others. Details on how we fit the emulators and demonstrate their accuracy can be found in Sections .2 and .3.

We are now in a position to build an inference algorithm using the site-specific stochastic emulators instead of the deterministic wavefront model . In doing so, we fix each emulator’s hyperparameters at their posterior mean. It is possible to undertake inferences allowing for uncertainty in the hyperparameters. We describe methods for doing this in Section .4 and also show that allowing for such uncertainty in our analysis has only a very minor affect on the posterior distribution and hence on our inferences.

Using the subscript to denote the use of the emulators, the (joint) posterior density for the unknown parameters in the statistical model, replacing (7), is

(12) |

where

has th element and, to allow for emulator uncertainty, we redefine , where . Sampling from the posterior distribution (12) can be achieved by using a Metropolis–Hastings scheme similar to that described in Section .2. Specifically, we use correlated normal random walks to propose updates on a log-scale in separate blocks for , and the parameters that govern spatial smoothness. The covariance structure of the innovations was chosen using a series of short trial runs.

## 6 Results

We found that running the MCMC scheme for 1M iterations after a burn-in of 100K iterations and thinning by taking every 100th iterate gave a sample of K (effectively uncorrelated) values from the posterior distribution. Figure 2 shows the marginal and joint posterior densities of the key parameters of the wavefront model (the diffusivity , the coastal velocity and the river-path velocity ) and also the parameters of the statistical model (the global error , and the Gaussian Process hyperparameters and ). The reduction in uncertainty from the prior to the posterior, which is considerable for most parameters, shows that the data have clearly been informative.

The posterior distributions of the wavefront model are in reasonable agreement with the range of values suggested by previous studies. For example, the mode of the marginal posterior for the diffusivity is around 20 kmyear, which corresponds to a global rate of spread of kmyear. The modal values for the amplitudes of the advective velocities, kmyear and kmyear, are rather lower than the values suggested in the archaeological literature which motivated their study [Zilhão (2001), Dolukhanov et al. (2005)]. Regarding the coastal speed , this could be explained by an accelerated spread only being required in the west Mediterranean, whereas our model applies it to the whole European coastline; an alternative model, allowing regional variations in coastal effects, may produce different results. The relatively small value for is perhaps more surprising, given the prevailing archaeological picture of relatively rapid spread of the LBK culture in the Danube and Rhine basins. It may be that the earliest dates associated with the LBK culture do not correspond particularly well to spread as a continuous wave and that an alternative model might also be preferable here. After taking account of the overall level of fit of our current model, the dates within this region are satisfactorily explained with a relatively weak advective enhancement of the background diffusive spread. It is also plausible that the slightly enhanced value of (compared with many of the studies cited in Section 1) results in the relatively low values for and .

The marginal posterior distribution for the global error parameter is centred on a value of around 575 years (see Figure 2). This value quantifies the mismatch between our mathematical model of the spread and the local variations present in the actual spread. Thus, our inference suggests that a simple large-scale wave of advance, while remaining a good model on the continental scale, should only be considered a good model on time scales of order 600 years (and thus length scales of order 600 km) or greater; on shorter time scales, significant local variations should be expected. Mismatch between model and data is also allowed in our spatial term, ; the marginal posterior distributions for its amplitude () and length scale () are relatively tightly peaked around values of order 30 years and , respectively. The latter value suggests that correlations between dates at nearby sites are effectively significant within distances of order 750 km. (Note that this is approximately the length scale on which suggests that the variations should be expected within the large scale model, so that the two results are consistent.) However, the amplitude is relatively small compared with , suggesting that “within region” variation is not particularly significant, in comparison with larger scale deviations from the model.

We can examine the operation of the spatial Gaussian process in more detail by looking at its distribution at various sites. Now

(13) |

where

and so posterior realizations from the spatial process can be obtained using the MCMC output , , by simulating from , . For example, the posterior spatial process at the UK site Monamore (W, N) has mean 21.7 years and standard deviation 33.0 years, and this small mean reflects both the good fit of the wavefront model at this site and at its neighboring sites; see Figure 1. In contrast, the French site Greifensee (W, N) has mean -726 years and standard deviation 453 years, and this large mean highlights that the data from this site is inconsistent with its neighbors.

### 6.1 Model fit

We use predictive simulations to assess the validity of the statistical model and thereby the underlying model of the propagating wavefront. The (joint) posterior predictive distribution of the arrival times at the radiocarbon sites can be determined in a similar way to that for the (posterior) spatial process. First we have

where and is the spatial Gaussian process evaluated at the sites. We can integrate over in closed form using (13) to obtain

and thereby determine the predictive arrival time distribution by simulating realizations from , . Figure 3 shows the radiocarbon dates for four randomly chosen sites together with their posterior predictive

distributions determined using this simulation technique. These sites are at Monamore (W, N), Vrsnik-Tarinci (E, N), Galabnik (E, N) and Greifensee (W, N). The C dates cluster very close to the mode of the posterior densities in three cases and the discrepant site (d) is the one discussed in the previous section whose data is not consistent with its neighboring sites. Very encouragingly, we found similar clustering of site dates around their predictive mode for around 90% of the 302 sites. We conclude that the wavefront model (or more properly, the emulators) provides an adequate description of the data.

Figure 4 illustrates how the wavefront moves across Europe via a map of the means of the posterior predictive distributions (Figure 3) generated by interpolating the mean arrival time at each site onto a regular mesh, using cubic interpolation with the Matlab routine griddata. In the same figure we also display the standard deviation of the predictive distributions interpolated in the same manner, highlighting areas where there is large uncertainty in the timing of the wavefront arrival from our model.

## 7 Discussion

Our main goal in this paper has been to develop generic statistical tools rather than to propose a more complex model for the spread of the Neolithic than, for example, Ackland et al. (2007), Davison et al. (2009) and Isern and Fort (2010). Our approach to inferring parameters of the Neolithic expansion is based on a direct application of the wavefront model. To alleviate the computational demands of the MCMC scheme, we constructed Gaussian process emulators for the arrival time of the wavefront at each of the 302 radiocarbon-dated sites in our data set, with the input parameters obtained from the wavefront model. This approach has been shown to perform well and is likely to do so for other spatio-temporal models that have complex space–time dependencies and are slow to simulate.

This paper represents the first attempt of statistical inference for quantitative models of the Neolithic dispersal based on modern statistical methods, providing an opportunity for a more detailed and reliable analysis of the data than ever before. We have improved upon earlier model fitting for the Neolithic expansion by avoiding the introduction of a maximum precision of the Early Neolithic C age determinations at 100–200 years, obtained as the minimum empirical standard deviation of C dates in large homogeneous data sets for key Neolithic sites [Dolukhanov and Shukurov (2004), Dolukhanov et al. (2005), Davison et al. (2009)]. In fact, results presented here remain unaffected when such a maximum precision is introduced. This was not the case with the earlier analyses, and the robustness of the inference procedure suggested here is encouraging.

We believe that we have demonstrated the usefulness of our method as a generic tool and that adopting a statistically rigorous approach constitutes a significant advance on previous “wave of advance” models of the Neolithic. The posterior distribution of the parameters of our wavefront model () provides clearly bounded constraints on the human processes underlying the spread of farming. Also, the posterior distribution of the parameters of our statistical model () allows a quantitative assessment of the utility of our wavefront model (including the limits of applicability of such a large-scale model to smaller regions). The latter point has important implications for the archaeological interpretation of the Neolithic transition in Europe and should contribute to the ongoing debate about the characterization of the spread as a continuous phenomenon. For example, Rowley-Conwy (2011) argues that the expansion of farming was characterized by sporadic events, punctuated by long periods without outward expansion; Rowley-Conwy includes the spread of the LBK and Impressed Ware cultures as rapid events of this type. In light of the difficulties encountered here in modeling these spreads, noted in the previous section, it would be of interest to consider alternative models, which are perhaps better able to model the expansion within the relevant regions.

In terms of future work, the wavefront model can be improved in several respects. Most importantly, a more detailed description of regional variations appears to be necessary. In particular, C data suggest that the accelerated spread along the coastlines is most pronounced in the western Mediterranean and around the Iberian Peninsula, but not on the North Sea and Black Sea coasts. The approach presented here can be readily generalized to include such a variation and to provide a useful estimate or an upper limit on the coastal propagation speed in Northern Europe. Also, our understanding of the Neolithic spread can be enriched by adding other parameters to our model, such as the starting time and position of the source of the Neolithic expansion, and even additional sources as suggested by Davison et al. (2009), based on the C evidence and corroborated by Motuzaite-Matuzeviciute, Hunt and Jones (2009) using archaeobotanical evidence. Another major improvement would be to use the published C dates directly for inference, rather than combining them at each site as in (1). Model fits to the original radiocarbon data without preliminary selection, unavoidably involving subjective and possibly poorly justified criteria, have never been attempted in the studies of the Neolithic [Hazelwood and Steele (2004), Steele (2010)].

A perennial question in the ongoing debate about the nature of the Neolithic is the mechanism of the spread of the agropastoral lifestyle. The wavefront model used here applies equally well to demic spread, involving migration of people, and to cultural diffusion, that is, the transition to farming due to indigenous adaptation and knowledge transfer [Zvelebil, Domaní-ska and Dennell (1998)]. However, cultural diffusion may be modeled more explicitly, in terms of multiple, mutually-interacting populations; see, for example, Aoki, Shida and Shigesada (1996), Ackland et al. (2007) and Patterson et al. (2010). We intend to explore the fit of some of these other models and to determine how plausible these different models are as explanations of the radiocarbon data.

## Appendix

### .1 Experimental design

In order to fit the emulators, we first need to obtain the training data by running the wavefront model at a suitable choice of -values. Computing resources available to us allowed a -point design. Although using a regular lattice design is appealingly simple, it is not particularly efficient. Instead, we adopt a more commonly used design for fitting Gaussian processes, namely, the Latin Hypercube Design (LHD) [McKay, Beckman and Conover (1979)]. Designs of this class are space-filling in that they distribute points within a hypercube in parameter space more efficiently than a lattice design. Our 200-point LHD was constructed using the maximin option in the Matlab routine lhsdesign. We set the lower bounds of the hypercube to be the origin and used the upper one percentiles of the prior distribution as its upper bounds. The routine then simulated 5000 LHDs and took the design which maximized the minimum distance between design points. We then used a sequential strategy in which we repeatedly ran the inference algorithm (described in the following section), used the results to determine a conservative estimate of a hypercube containing all points in the MCMC output (and therefore plausibly containing all of the posterior density), and then generated another LHD. The final LHD used for inferences on in this paper is contained within the hypercube .

### .2 Fitting the emulator

At site , this essentially requires the determination of the posterior distribution of the emulator’s hyperparameters using the training data from the wavefront model evaluated at the emulator’s design points. We adopt fairly diffuse independent log-normal priors for the GP hyperparameters, with and . Then their (joint) posterior density is

where

For notational simplicity, we have dropped the explicit dependence of the training data , the mean and the covariance matrix on the training points . It is fairly straightforward to implement an MCMC scheme to simulate realizations from this posterior distribution. We used a Metropolis–Hastings scheme in which proposals for each hyperparameter were made on a log-scale using (univariate) normal random walks. A suitable magnitude of the innovations (their variance) was adopted after a series of short trial runs. For each site-specific emulator, the MCMC scheme was run for 500K iterations after a burn-in of 50K iterations and the output thinned by taking every 50th iterate, giving a total of 10K (almost uncorrelated) iterates to be used for inference. Figure 5 shows the marginal posterior densities of the emulator hyperparameters for site 1, the Achilleion

site, located at latitude N, longitude E. Note that these posterior densities have relatively small variability and are robust to the parameter choice for the prior of the hyperparameters, and both these features are typical of all the sites in our data set.

### .3 Emulator validation

The accuracy of our fitted emulators as an approximation to the wavefront model can be assessed by using a variety of graphical and numerical methods. Our assessment was based on the fit at a set of parameter values that were not used to fit the emulator: we used a -point Latin hypercube design . First, for each site , we obtained the arrival times from the propagating front model. We then obtained realizations of the emulator at these design points which properly accounted for the uncertainty in the hyperparameters by using (9) and the realizations from the MCMC fit of the emulator. Plots of the emulator’s mean and 95% credible region showed that the emulators provided a reasonably accurate fit to the wavefront output. We note that these plots were very similar to those determined by fixing the hyperparameters at their posterior means.

We also assessed the fit of our emulators by using a numerical statistic taken from Bastos and O’Hagan (2009) which compares the emulator output with that from the wavefront model allowing for the site-specific accuracy and correlation between design points. This statistic is the site-specific Mahalanobis distance MD, , where

which has a scaled -distribution in the case of the GP emulator with fixed hyperparameters, with . In our calculations we have fixed the hyperparameters at their posterior means. Figure 6 shows the Mahalanobis distance at each

site, together with the upper point of its distribution, and the spatial distribution of the Mahalanobis distance. The figure confirms that the emulators provide a reasonable fit throughout the design space and that there is no obvious spatial dependence in the fit of the emulators. We note that these diagnostics gave similar conclusions for all values in the posterior sample of and for different LHDs without any systematic site-specific biases.

Finally, we have assessed the distributional fit of the emulators by calculating the probability integral transform (PIT) values [Gneiting, Balabdaoui and Raftery (2007)] for each site at each of the 100 points in the design. Again, we fix the hyperparameters at their posterior means. A well fitting emulator should produce values from the standard uniform distribution. A histogram of the PIT values over all sites is shown in Figure 7. The figure also shows a summary of how closely the PIT values from each site resemble the standard uniform distribution. Here we use a Pearson statistic after binning the data into ten equal sized bins and so also show the 95% point of the distribution. This figure confirms that the emulator fits are satisfactory.

### .4 Allowing for hyperparameter uncertainty

In Section 5 we describe how to obtain realizations from the posterior distribution (12) when fixing the emulator hyperparameters at their posterior mean. Here we extend the method to take account of hyperparameter uncertainty by using their realizations obtained when fitting the emulators.

We first describe an MCMC scheme in which the emulators themselves are marginalized over hyperparameter uncertainty. In this case, the marginal distribution of can be well approximated by an equally weighted normal mixture, namely,

where and are given by (10) and (11). Therefore, the (emulator) likelihood in (12) can be written as

where is the density of an -dimensional normal distribution with mean and covariance matrix , has th element and, to allow for the emulator uncertainty, we redefine , where . Sampling from the posterior distribution when using this (emulator) likelihood can be achieved by using the same methods as described in Section 5. However, a major drawback of using this scheme is the time taken to calculate the (emulator) likelihood at each iteration. Therefore, we prefer to use an MCMC scheme which uses a joint update for the model parameters and the hyperparameters and has (12) as a marginal distribution. Specifically, in the joint update, the value proposed for the hyperparameters is sampled from the posterior realizations described above. Although this scheme has a much larger state space, each MCMC iteration is considerably faster than in the marginal scheme and we have found it to have good convergence properties.

We ran the MCMC scheme for 1M iterations (after a burn-in of 100K iterations), and then thinned the output taking every 100th iterate to leave an effectively uncorrelated sample of 10K draws from the posterior densities. Figure 8 shows the marginal posterior distributions allowing for emulator hyperparameter uncertainty and for fixed hyperparameters. We note that it is difficult to distinguish between these cases. Similar conclusions were found (when comparing fixed with uncertain hyperparameters) for the site-specific posterior distribution of the spatial process and for the predictive distributions .

## Acknowledgments

The authors would like to thank an Editor and three anonymous referees for their constructive comments on an earlier draft.

Data \slink[doi]10.1214/12-AOAS579SUPP \slink[url]http://lib.stat.cmu.edu/aoas/579/supplement.zip \sdatatype.zip \sdescriptionThe zip file contains the raw data (xls file) and a description of the raw data (pdf file).

## References

- Ackland et al. (2007) {barticle}[pbm] \bauthor\bsnmAckland, \bfnmGraeme J.\binitsG. J., \bauthor\bsnmSignitzer, \bfnmMarkus\binitsM., \bauthor\bsnmStratford, \bfnmKevin\binitsK. \AND\bauthor\bsnmCohen, \bfnmMorrel H.\binitsM. H. (\byear2007). \btitleCultural hitchhiking on the wave of advance of beneficial technologies. \bjournalProc. Natl. Acad. Sci. USA \bvolume104 \bpages8714–8719. \biddoi=10.1073/pnas.0702469104, issn=0027-8424, pii=0702469104, pmcid=1885568, pmid=17517663 \bptokimsref \endbibitem
- Aitken (1990) {bbook}[author] \bauthor\bsnmAitken, \bfnmM.\binitsM. (\byear1990). \btitleScience-Based Dating in Archaeology. \bpublisherLongman, \blocationHarlow. \bptokimsref \endbibitem
- Ammerman and Cavalli-Sforza (1971) {barticle}[author] \bauthor\bsnmAmmerman, \bfnmA. J.\binitsA. J. \AND\bauthor\bsnmCavalli-Sforza, \bfnmL. L.\binitsL. L. (\byear1971). \btitleMeasuring the rate of spread of early farming in Europe. \bjournalMan \bvolume6 \bpages674–688. \bptokimsref \endbibitem
- Ammerman and Cavalli-Sforza (1973) {bincollection}[author] \bauthor\bsnmAmmerman, \bfnmA. J.\binitsA. J. \AND\bauthor\bsnmCavalli-Sforza, \bfnmL. L.\binitsL. L. (\byear1973). \btitleA population model for the diffusion of early farming in Europe. In \bbooktitleThe Explanation of Culture Change (\beditor\bfnmColin\binitsC. \bsnmRenfrew, ed.) \bpages343–357. \bpublisherDuckworth, \blocationLondon. \bptokimsref \endbibitem
- Aoki, Shida and Shigesada (1996) {barticle}[author] \bauthor\bsnmAoki, \bfnmKenichi\binitsK., \bauthor\bsnmShida, \bfnmMitsuo\binitsM. \AND\bauthor\bsnmShigesada, \bfnmNanako\binitsN. (\byear1996). \btitleTravelling wave solutions for the spread of farmers into a region occupied by hunter-gatherers. \bjournalTheor. Popul. Biol. \bvolume50 \bpages1–17. \bptokimsref \endbibitem
- Baggaley et al. (2009) {barticle}[author] \bauthor\bsnmBaggaley, \bfnmAndrew W.\binitsA. W., \bauthor\bsnmBarenghi, \bfnmCarlo F.\binitsC. F., \bauthor\bsnmShukurov, \bfnmAnvar\binitsA. \AND\bauthor\bsnmSubramanian, \bfnmKandaswamy\binitsK. (\byear2009). \btitleReconnecting flux-rope dynamo. \bjournalPhys. Rev. E \bvolume80 \bpages055301. \bptokimsref \endbibitem
- Baggaley et al. (2013) {bmisc}[author] \bauthor\bsnmBaggaley, \bfnmAndrew W.\binitsA. W., \bauthor\bsnmBoys, \bfnmRichard J.\binitsR. J., \bauthor\bsnmGolightly, \bfnmAndrew\binitsA. \AND\bauthor\bsnmShukurov, \bfnmAnvar\binitsA. (\byear2013). \bhowpublishedSupplement to “Inference for population dynamics in the Neolithic period.” DOI:\doiurl10.1214/12-AOAS579SUPP. \bptokimsref \endbibitem
- Bastos and O’Hagan (2009) {barticle}[mr] \bauthor\bsnmBastos, \bfnmLeonardo S.\binitsL. S. \AND\bauthor\bsnmO’Hagan, \bfnmAnthony\binitsA. (\byear2009). \btitleDiagnostics for Gaussian process emulators. \bjournalTechnometrics \bvolume51 \bpages425–438. \biddoi=10.1198/TECH.2009.08019, issn=0040-1706, mr=2756478 \bptokimsref \endbibitem
- Bayarri et al. (2007) {barticle}[mr] \bauthor\bsnmBayarri, \bfnmMaria J.\binitsM. J., \bauthor\bsnmBerger, \bfnmJames O.\binitsJ. O., \bauthor\bsnmPaulo, \bfnmRui\binitsR., \bauthor\bsnmSacks, \bfnmJerry\binitsJ., \bauthor\bsnmCafeo, \bfnmJohn A.\binitsJ. A., \bauthor\bsnmCavendish, \bfnmJames\binitsJ., \bauthor\bsnmLin, \bfnmChin-Hsu\binitsC.-H. \AND\bauthor\bsnmTu, \bfnmJian\binitsJ. (\byear2007). \btitleA framework for validation of computer models. \bjournalTechnometrics \bvolume49 \bpages138–154. \biddoi=10.1198/004017007000000092, issn=0040-1706, mr=2380530 \bptokimsref \endbibitem
- Birdsell (1957) {barticle}[author] \bauthor\bsnmBirdsell, \bfnmJ. B.\binitsJ. B. (\byear1957). \btitleSome population problems involving Pleistocene man. \bjournalCold Spring Harbor Symposium on Quantitative Biology \bvolume22 \bpages47–69. \bptokimsref \endbibitem
- Blackwell and Buck (2008) {barticle}[mr] \bauthor\bsnmBlackwell, \bfnmP. G.\binitsP. G. \AND\bauthor\bsnmBuck, \bfnmC. E.\binitsC. E. (\byear2008). \btitleEstimating radiocarbon calibration curves. \bjournalBayesian Anal. \bvolume3 \bpages225–248. \bidissn=1936-0975, mr=2407423 \bptokimsref \endbibitem
- Bocquet-Appel et al. (2009) {barticle}[author] \bauthor\bsnmBocquet-Appel, \bfnmJean-Pierré\binitsJ.-P., \bauthor\bsnmNaji, \bfnmStephan\binitsS., \bauthor\bsnmLinden, \bfnmMarc Vander\binitsM. V. \AND\bauthor\bsnmKozlowski, \bfnmJanusz K.\binitsJ. K. (\byear2009). \btitleDetection of diffusion and contact zones of early farming in Europe from the space–time distribution of 14C dates. \bjournalJ. Archeo. Sci. \bvolume36 \bpages807–820. \bptokimsref \endbibitem
- Buck et al. (1991) {barticle}[author] \bauthor\bsnmBuck, \bfnmC. E.\binitsC. E., \bauthor\bsnmKenworthy, \bfnmJ. B.\binitsJ. B., \bauthor\bsnmLitton, \bfnmC. D.\binitsC. D. \AND\bauthor\bsnmSmith, \bfnmA. F. M.\binitsA. F. M. (\byear1991). \btitleCombining archaeological and radiocarbon information: A Bayesian approach to calibration. \bjournalAntiquity \bvolume65 \bpages808–821. \bptokimsref \endbibitem
- Burger and Thomas (2011) {bincollection}[author] \bauthor\bsnmBurger, \bfnmJoachim\binitsJ. \AND\bauthor\bsnmThomas, \bfnmMark G.\binitsM. G. (\byear2011). \btitleThe palaeopopulationgenetics of humans, cattle and diarying in neolithic Europe. In \bbooktitleHuman Bioarchaeology of the Transition to Agriculture (\beditor\bfnmRon\binitsR. \bsnmPinhasi \AND\beditor\bfnmJay T.\binitsJ. T. \bsnmStock, eds.) \bpages1029–1058. \bpublisherWiley, \blocationChichester. \bptokimsref \endbibitem
- Davison et al. (2006) {barticle}[author] \bauthor\bsnmDavison, \bfnmK.\binitsK., \bauthor\bsnmDolukhanov, \bfnmP.\binitsP., \bauthor\bsnmSarson, \bfnmG. R.\binitsG. R. \AND\bauthor\bsnmShukurov, \bfnmA.\binitsA. (\byear2006). \btitleThe role of waterways in the spread of the Neolithic. \bjournalJ. Archeo. Sci. \bvolume33 \bpages641–652. \bptokimsref \endbibitem
- Davison et al. (2009) {barticle}[author] \bauthor\bsnmDavison, \bfnmK.\binitsK., \bauthor\bsnmDolukhanov, \bfnmP. M.\binitsP. M., \bauthor\bsnmSarson, \bfnmG. R.\binitsG. R., \bauthor\bsnmShukurov, \bfnmA.\binitsA. \AND\bauthor\bsnmZaitseva, \bfnmG. I.\binitsG. I. (\byear2009). \btitleMultiple sources of the European Neolithic: Mathematical modelling constrained by radiocarbon dates. \bjournalQuaternary International \bvolume203 \bpages10–18. \bptokimsref \endbibitem
- Dolukhanov and Shukurov (2004) {barticle}[author] \bauthor\bsnmDolukhanov, \bfnmPavel\binitsP. \AND\bauthor\bsnmShukurov, \bfnmAnvar\binitsA. (\byear2004). \btitleModelling the Neolithic dispersal in Northern Eurasia. \bjournalDocumenta Praehistorica \bvolume31 \bpages35–47. \bptokimsref \endbibitem
- Dolukhanov et al. (2005) {barticle}[author] \bauthor\bsnmDolukhanov, \bfnmPavel\binitsP., \bauthor\bsnmShukurov, \bfnmAnvar\binitsA., \bauthor\bsnmGronenborn, \bfnmDetlef\binitsD., \bauthor\bsnmSokoloff, \bfnmDmitry\binitsD., \bauthor\bsnmTimofeev, \bfnmVladimir\binitsV. \AND\bauthor\bsnmZaitseva, \bfnmGanna\binitsG. (\byear2005). \btitleThe chronology of Neolithic dispersal in Central and Eastern Europe. \bjournalJ. Archeo. Sci. \bvolume32 \bpages1441–1458. \bptokimsref \endbibitem
- Fisher (1937) {barticle}[author] \bauthor\bsnmFisher, \bfnmR. A.\binitsR. A. (\byear1937). \btitleThe wave of advance of advantageous genes. \bjournalAnn. Eugenics \bvolume7 \bpages353–369. \bptokimsref \endbibitem
- Fort and Pujol (2008) {barticle}[author] \bauthor\bsnmFort, \bfnmJoaquim\binitsJ. \AND\bauthor\bsnmPujol, \bfnmPujol Toni\binitsP. T. (\byear2008). \btitleProgress in front propagation research. \bjournalRep. Progr. Phys. \bvolume71 \bpages086001 (41 pp.). \bptokimsref \endbibitem
- Geophysical Data System (2011) {bmisc}[author] \borganizationGeophysical Data System (\byear2011). \bhowpublishedAvailable at http://www.ngdc.noaa.gov/mgg/geodas/ geodas.html. \bptokimsref \endbibitem
- Gkiasta et al. (2003) {barticle}[author] \bauthor\bsnmGkiasta, \bfnmMarina\binitsM., \bauthor\bsnmRussell, \bfnmThembi\binitsT., \bauthor\bsnmShennan, \bfnmStephen\binitsS. \AND\bauthor\bsnmSteele, \bfnmJames\binitsJ. (\byear2003). \btitleNeolithic transition in Europe: The radiocarbon record revisited. \bjournalAntiquity \bvolume77 \bpages45–62. \bptokimsref \endbibitem
- Gneiting, Balabdaoui and Raftery (2007) {barticle}[mr] \bauthor\bsnmGneiting, \bfnmTilmann\binitsT., \bauthor\bsnmBalabdaoui, \bfnmFadoua\binitsF. \AND\bauthor\bsnmRaftery, \bfnmAdrian E.\binitsA. E. (\byear2007). \btitleProbabilistic forecasts, calibration and sharpness. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume69 \bpages243–268. \biddoi=10.1111/j.1467-9868.2007.00587.x, issn=1369-7412, mr=2325275 \bptokimsref \endbibitem
- Hazelwood and Steele (2004) {barticle}[author] \bauthor\bsnmHazelwood, \bfnmLee\binitsL. \AND\bauthor\bsnmSteele, \bfnmJames\binitsJ. (\byear2004). \btitleSpatial dynamics of human dispersals. Constraints on modelling and archaeological validation. \bjournalJ. Archeo. Sci. \bvolume31 \bpages669–679. \bptokimsref \endbibitem
- Henderson et al. (2009) {barticle}[mr] \bauthor\bsnmHenderson, \bfnmDaniel A.\binitsD. A., \bauthor\bsnmBoys, \bfnmRichard J.\binitsR. J., \bauthor\bsnmKrishnan, \bfnmKim J.\binitsK. J., \bauthor\bsnmLawless, \bfnmConor\binitsC. \AND\bauthor\bsnmWilkinson, \bfnmDarren J.\binitsD. J. (\byear2009). \btitleBayesian emulation and calibration of a stochastic computer model of mitochondrial DNA deletions in substantia nigra neurons. \bjournalJ. Amer. Statist. Assoc. \bvolume104 \bpages76–87. \biddoi=10.1198/jasa.2009.0005, issn=0162-1459, mr=2663034 \bptokimsref \endbibitem
- Isern and Fort (2010) {barticle}[author] \bauthor\bsnmIsern, \bfnmN.\binitsN. \AND\bauthor\bsnmFort, \bfnmJ.\binitsJ. (\byear2010). \btitleAnisotropic dispersion, space competition and the slowdown of the Neolithic transition. \bjournalNew J. Phys. \bvolume12 \bpages123002. \bptokimsref \endbibitem
- Kennedy and O’Hagan (2001) {barticle}[mr] \bauthor\bsnmKennedy, \bfnmMarc C.\binitsM. C. \AND\bauthor\bsnmO’Hagan, \bfnmAnthony\binitsA. (\byear2001). \btitleBayesian calibration of computer models. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume63 \bpages425–464. \biddoi=10.1111/1467-9868.00294, issn=1369-7412, mr=1858398 \bptokimsref \endbibitem
- Kolmogorov, Petrovskii and Piskunov (1937) {barticle}[author] \bauthor\bsnmKolmogorov, \bfnmA.\binitsA., \bauthor\bsnmPetrovskii, \bfnmI.\binitsI. \AND\bauthor\bsnmPiskunov, \bfnmN.\binitsN. (\byear1937). \btitleA study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem. \bjournalByul. Moskovskogo Gos. Univ. \bvolume1 \bpages1–25. \bptokimsref \endbibitem
- McKay, Beckman and Conover (1979) {barticle}[mr] \bauthor\bsnmMcKay, \bfnmM. D.\binitsM. D., \bauthor\bsnmBeckman, \bfnmR. J.\binitsR. J. \AND\bauthor\bsnmConover, \bfnmW. J.\binitsW. J. (\byear1979). \btitleA comparison of three methods for selecting values of input variables in the analysis of output from a computer code. \bjournalTechnometrics \bvolume21 \bpages239–245. \biddoi=10.2307/1268522, issn=0040-1706, mr=0533252 \bptokimsref \endbibitem
- Motuzaite-Matuzeviciute, Hunt and Jones (2009) {bincollection}[author] \bauthor\bsnmMotuzaite-Matuzeviciute, \bfnmG.\binitsG., \bauthor\bsnmHunt, \bfnmH. V.\binitsH. V. \AND\bauthor\bsnmJones, \bfnmM. K.\binitsM. K. (\byear2009). \btitleMultiple sources for Neolithic European agriculture: Geographical origins of early domesticates in Moldova and Ukraine. In \bbooktitleThe East European Plain on the Eve of Agriculture (\beditor\bfnmPavel M.\binitsP. M. \bsnmDolukhanov, \beditor\bfnmGraeme R.\binitsG. R. \bsnmSarson \AND\beditor\bfnmAnvar\binitsA. \bsnmShukurov, eds.). \bseriesBritish Archaeological Reports, International Series \bvolume1964 \bpages53–64. \bpublisherArchaeopress, \blocationOxford. \bptokimsref \endbibitem
- Pape (2004) {bmisc}[author] \bauthor\bsnmPape, \bfnmD.\binitsD. (\byear2004). \bhowpublishedWorld DataBank II. Available at http://www.evl.uic.edu/pape/data/ WDB/. \bptokimsref \endbibitem
- Patterson et al. (2010) {barticle}[author] \bauthor\bsnmPatterson, \bfnmM. A.\binitsM. A., \bauthor\bsnmSarson, \bfnmG. R.\binitsG. R., \bauthor\bsnmSarson, \bfnmH. C.\binitsH. C. \AND\bauthor\bsnmShukurov, \bfnmA.\binitsA. (\byear2010). \btitleModelling the Neolithic transition in a heterogeneous environment. \bjournalJ. Archeo. Sci. \bvolume37 \bpages2929–2937. \bptokimsref \endbibitem
- Pettitt et al. (2003) {barticle}[author] \bauthor\bsnmPettitt, \bfnmP. B.\binitsP. B., \bauthor\bsnmDavies, \bfnmW.\binitsW., \bauthor\bsnmGamble, \bfnmC. S.\binitsC. S. \AND\bauthor\bsnmRichards, \bfnmM. B.\binitsM. B. (\byear2003). \btitlePalaeolithic radiocarbon chronology: Quantifying our confidence beyond two half-lives. \bjournalJ. Archaeo. Sci. \bvolume30 \bpages1685–1693. \bptokimsref \endbibitem
- Rasmussen and Williams (2006) {bbook}[mr] \bauthor\bsnmRasmussen, \bfnmCarl Edward\binitsC. E. \AND\bauthor\bsnmWilliams, \bfnmChristopher K. I.\binitsC. K. I. (\byear2006). \btitleGaussian Processes for Machine Learning. \bpublisherMIT Press, \blocationCambridge, MA. \bidmr=2514435 \bptokimsref \endbibitem
- Reimer et al. (2004) {barticle}[author] \bauthor\bsnmReimer, \bfnmP. J.\binitsP. J., \bauthor\bsnmBaillie, \bfnmM.\binitsM., \bauthor\bsnmBard, \bfnmE.\binitsE., \bauthor\bsnmBayliss, \bfnmA.\binitsA., \bauthor\bsnmBeck, \bfnmJ.\binitsJ., \bauthor\bsnmBertrand, \bfnmC.\binitsC., \bauthor\bsnmBlackwell, \bfnmP.\binitsP., \bauthor\bsnmBuck, \bfnmC.\binitsC., \bauthor\bsnmBurr, \bfnmG.\binitsG., \bauthor\bsnmCutler, \bfnmK.\binitsK., \bauthor\bsnmDamon, \bfnmP.\binitsP., \bauthor\bsnmEdwards, \bfnmR.\binitsR., \bauthor\bsnmFairbanks, \bfnmR.\binitsR., \bauthor\bsnmFriedrich, \bfnmM.\binitsM., \bauthor\bsnmGuilderson, \bfnmT.\binitsT., \bauthor\bsnmHughen, \bfnmK.\binitsK., \bauthor\bsnmKromer, \bfnmB.\binitsB., \bauthor\bsnmMcCormac, \bfnmF.\binitsF., \bauthor\bsnmManning, \bfnmS.\binitsS., \bauthor\bsnmBronk Ramsey, \bfnmC.\binitsC., \bauthor\bsnmReimer, \bfnmR.\binitsR., \bauthor\bsnmRemmele, \bfnmS.\binitsS., \bauthor\bsnmSouthon, \bfnmJ.\binitsJ., \bauthor\bsnmStuiver, \bfnmM.\binitsM., \bauthor\bsnmTalamo, \bfnmS.\binitsS., \bauthor\bsnmTaylor, \bfnmF.\binitsF., \bauthor\bparticlevan der \bsnmPlicht, \bfnmJ.\binitsJ. \AND\bauthor\bsnmWeyhenmeyer, \bfnmC.\binitsC. (\byear2004). \btitleIntCal04 terrestrial radiocarbon age calibration. \bjournalRadiocarbon \bvolume46 \bpages1029–1058. \bptokimsref \endbibitem
- Rowley-Conwy (2011) {barticle}[author] \bauthor\bsnmRowley-Conwy, \bfnmPeter\binitsP. (\byear2011). \btitleWestward ho! The spread of agriculture from Central Europe to the Atlantic. \bjournalCurrent Anthropology \bvolume52 \bpagesS431–S451. \bptokimsref \endbibitem
- Santner, Williams and Notz (2003) {bbook}[mr] \bauthor\bsnmSantner, \bfnmThomas J.\binitsT. J., \bauthor\bsnmWilliams, \bfnmBrian J.\binitsB. J. \AND\bauthor\bsnmNotz, \bfnmWilliam I.\binitsW. I. (\byear2003). \btitleThe Design and Analysis of Computer Experiments. \bpublisherSpringer, \blocationNew York. \bidmr=2160708 \bptokimsref \endbibitem
- Scott, Cook and Naysmith (2007) {barticle}[author] \bauthor\bsnmScott, \bfnmE. Marian\binitsE. M., \bauthor\bsnmCook, \bfnmGordon T.\binitsG. T. \AND\bauthor\bsnmNaysmith, \bfnmPhilip\binitsP. (\byear2007). \btitleError and uncertainty in radiocarbon measurements. \bjournalRadiocarbon \bvolume49 \bpages427–440. \bptokimsref \endbibitem
- Steele (2010) {barticle}[author] \bauthor\bsnmSteele, \bfnmJames\binitsJ. (\byear2010). \btitleRadiocarbon dates as data: Quantitative strategies for estimating colonization front speeds and event densities. \bjournalJ. Archeo. Sci. \bvolume37 \bpages2017–2030. \bptokimsref \endbibitem
- Steele, Adams and Sluckin (1998) {barticle}[author] \bauthor\bsnmSteele, \bfnmJames\binitsJ., \bauthor\bsnmAdams, \bfnmJonathan\binitsJ. \AND\bauthor\bsnmSluckin, \bfnmTim\binitsT. (\byear1998). \btitleModelling Paleoindian dispersals. \bjournalWorld Archaeology \bvolume30 \bpages286–305. \bptokimsref \endbibitem
- Whittle (1996) {bbook}[author] \bauthor\bsnmWhittle, \bfnmAlasdair\binitsA. (\byear1996). \btitleEurope in the Neolithic: The Creation of New Worlds. \bpublisherCambridge Univ. Press, \blocationCambridge. \bptokimsref \endbibitem
- Zilhão (2001) {barticle}[author] \bauthor\bsnmZilhão, \bfnmJoao\binitsJ. (\byear2001). \btitleRadiocarbon evidence for maritime pioneer colonization at the origins of farming in west Mediterranean Europe. \bjournalProc. Natl. Acad. Sci. USA \bvolume98 \bpages14180–14185. \bptokimsref \endbibitem
- Zvelebil, Domaní-ska and Dennell (1998) {bbook}[author] \bauthor\bsnmZvelebil, \bfnmM.\binitsM., \bauthor\bsnmDomaní-ska, \bfnmL.\binitsL. \AND\bauthor\bsnmDennell, \bfnmR.\binitsR. (\byear1998). \btitleHarvesting the Sea, Farming the Forest: The Emergence of Neolithic Societies in the Baltic Region. \bpublisherSheffield Academic Press, \blocationSheffield. \bptokimsref \endbibitem