Therefore, it is important to ensure the date range of the data and model are appropriate for each other, and to exclude dates from the dataset that do not reasonably fall within the modelled range. We achieve this with our real datasets by only including a date if more than 50% of its probability falls within the modelled date range-i.e. it is more probable that its true date is internal than external. Similarly, we achieve this with our extremely small toy dataset (N = 6) by constraining the modelled date range to exclude the negligible tails outside the calibrated dates.
7. Search algorithm for parameters
The CPL model is a PMF such that the probability outside the date range equals 0, and the total probability within the date range equals 1. The exact shape of this PMF is defined by the (x, y) coordinates of the hinge points. Therefore, there are various constraints on parameters required to define such a curve. For example, if we consider a 2-CPL model, only the middle hinge has a free x-coordinate parameter, since the start and end date are already specified by the date range. Of the three y-coordinates (left, middle, right hinges), only two are free parameters, since the total probability must equal 1. Therefore, a 2-CPL model has three free parameters (one x-coordinate and two y-coordinates) and an n-phase CPL model has 2n?1 free parameters.
We perform the search for the ML parameters (given a 14 C dataset and calibration curve) using the differential evolution optimization algorithm DEoptimR . A naive approach to this search would propose a set of values for all parameters in an iteration simultaneously, and reject the set if it does not satisfy the above constraints. However, this approach would result in the rejection of many parameter sets. Instead, our objective function considers the parameters in order, such that the next parameter is searched for in a reduced parameter space, conditional on the previous parameters. We achieve this by adapting the ‘stick breaking’ Dirichlet process to apply in two dimensions by sampling stick breaks on the x-axis using the beta distribution and y-coordinates using the gamma distribution. At each hinge, the length of the stick is constrained by calculating the total area so far between the first and previous hinge.
Having constructed a likelihood function that calculates the relative likelihood of any parameter combination, 
it can be used as the objective function in a parameter search to find the ML parameter estimates. However, we also use the likelihood function in a ework to estimate credible intervals of our parameter estimates. We achieve this using the Metropolis–Hastings algorithm using a single chain of 100 000 iterations, discarding the first 2000 for burn-in, and thinning to every fifth iteration. The resulting joint posterior distribution can then be graphically represented in several ways, such as histograms of the marginal distributions (figure 6) or directly plotting the joint parameter estimates on a two-dimensional plot (figure 7).
9. Goodness-of-fit test
Once the best CPL model has been selected, its parameters found and the likelihood calculated, we generate 1000 simulated 14 C datasets under this CPL model by ‘uncalibrating’ calendar dates randomly sampled under the model, taking care to ensure sample sizes exactly match the number of phases in the observed dataset. We then calculate the proportion of each calibrated simulated dataset outside the 95% CI, giving a distribution of summary statistics under our best CPL model. The p-value is then calculated as the proportion of these simulated summary statistics that are smaller or equal to the observed summary statistic. Conceptually, this is similar to the method of calculating p-values under existing simulation methods for testing a null model [12,25–33].
