2.1 Activity pulses and Monte Carlo KDE–minima algorithm

Let $\{ {p_i}\left( t \right)\} _{i = 1}^n$ denote calibrated posterior distributions for $n$ radiocarbon determinations, where $t$ is calendar age. The goal is to partition the time axis into contiguous intervals within which the underlying event intensity can reasonably be treated as approximately homogeneous, and between which there is evidence for a lull in activity. These intervals are termed activity pulses. They are explicitly scale-dependent features of a smoothed intensity estimate at a chosen bandwidth $h$ and grid resolution, not direct estimates of a unique demographic trajectory (Silverman Reference Silverman2018; Scott Reference Scott2015).

The Monte Carlo KDE–minima procedure is:

1. 1. For each replicate $r = 1, \ldots, R$ , draw a single calendar age $t_i^{\left( r \right)}$ from each posterior ${p_i}\left( t \right)$ . Across replicates, these draws approximate integrating over the full posteriors (Efron and Tibshirani Reference Efron and Tibshirani1994; Bronk Ramsey Reference Bronk Ramsey2017).

2. 2. For each replicate, construct a KDE

   ${{\hat f}^{\left( r \right)}}\left( t \right) = {1 \over {nh}}\sum\limits_{i = 1}^n K \left( {{{t - t_i^{\left( r \right)}} \over h}} \right)$
   on a fixed calendar grid, where $K$ is the standard normal density and $h$ the bandwidth. The grid spacing is fine enough that minima can be located at a resolution of a few years relative to $h$ (Silverman Reference Silverman2018; Scott Reference Scott2015).

3. 3. At each grid point ${t_j}$ , compute the pointwise median $\tilde f\left( {{t_j}} \right)$ and the 2.5–97.5% Monte Carlo envelope $\left[ {{{\tilde f}_{0.025}}\left( {{t_j}} \right),{{\tilde f}_{0.975}}\left( {{t_j}} \right)} \right]$ over replicates.

4. 4. Identify local minima of $\tilde f\left( t \right)$ and of the lower envelope ${\tilde f_{0.025}}\left( t \right)$ using the sign pattern of first differences along the grid.

5. 5. Mark a grid index ${j^{\rm{*}}}$ as a persistent minimum if it is a local minimum of the median KDE and there exists a minimum of the lower envelope within one grid step of ${j^{\rm{*}}}$ .

6. 6. Let accepted minima occur at calendar positions ${b_1}{\rm{ }} \lt \ldots{\rm{ }} \lt {b_{K - 1}}$ . Define $K$ activity pulses as the intervals $\left( { - \infty, {b_1}} \right],\left( {{b_1},{b_2}} \right], \ldots, \left( {{b_{K - 1}}, + \infty } \right)$ , truncated to the calendar range spanned by the data.

The requirement that minima be present in both the median and the lower envelope is intentionally conservative. A minimum visible in $\tilde f\left( t \right)$ but absent from ${\tilde f_{0.025}}\left( t \right)$ is treated as a likely artefact of sampling variability or calibration-induced shoulders. The single-grid-step neighbourhood allows for minor horizontal shifts in minima between the median and envelope. In spirit, the procedure is analogous to multi-scale curve exploration methods such as SiZer, which emphasise persistent structures over those that disappear under small perturbations (Chaudhuri and Marron Reference Chaudhuri and Marron1999). The examples in Figure 1 show how the algorithm behaves in an “easy” case with widely separated pulses and a “hard” case with fewer dates and substantially overlapping episodes.

In all simulations and empirical analyses, calendar grids are restricted to the window of interest, and pulses abutting the edges of this window must be interpreted cautiously because minima near boundaries are less stable. In practice, analysts should examine whether edge pulses persist when the analysis window is expanded or when $h$ is varied within a plausible range (Silverman Reference Silverman2018; Scott Reference Scott2015). The bandwidth range used throughout this study, $h = 80$ –120 yr, should be read as a decision-relevant multi-century smoothing regime rather than as an attempt to estimate an optimal bandwidth. The substantive target is the recovery of broad, multi-century pulses, not the finest-scale local structure visible in an SPD or KDE curve. Given the 400-year pulse widths used in the simulations, values in this range are large enough to suppress small shoulders and calibration-induced irregularities but still small enough not to force adjacent 400-year pulses to merge a priori. Formal bandwidth selectors exist (Sheather Reference Sheather2004; Silverman Reference Silverman2018; Scott Reference Scott2015), but for radiocarbon data with irregular calibration structure and posterior uncertainty, sensitivity over a plausible range is more informative than any single plug-in estimate. Bandwidth sensitivity for the KRT application is examined in the empirical illustration (Section 3.3).

2.2 Level A simulations: calendar-time pulses

Level A simulations operate directly in calendar time (BCE), isolating the effects of sample size, episode spacing and bandwidth without calibration. The underlying intensity $\lambda \left( t \right)$ is piecewise constant with three pulses and a uniform background. The configuration is summarised in Table 1. In the baseline design, three pulses of width 400 yr are centred at $ - 3500$ , $ - 2700$ and $ - 1900$ cal BCE, that is

$\left[ {{s_1},{e_1}\left] = \right[ - 3700, - 3300\left] {,{\rm{\;\;}}} \right[{s_2},{e_2}\left] = \right[ - 2900, - 2500\left] {,{\rm{\;\;}}} \right[{s_3},{e_3}\left] = \right[ - 2100, - 1700} \right],$

with intensities ${\lambda _1} = 1$ , ${\lambda _2} = 1.5$ , ${\lambda _3} = 2$ and background intensity ${\lambda _{{\rm{bg}}}} = 0.1$ . The simulation window spans $ - 4000$ to $ - 1200$ cal BCE. Under these settings, the three pulses account for approximately 92% of expected events over the window. The low background component is included to prevent the simulations from collapsing into an unrealistically clean pulse-only system and to approximate the low-rate dated events that would often exist outside major activity clusters. This stylised, piecewise-constant configuration with equal-width pulses is not intended as a literal model of any particular archaeological sequence; it provides a controlled setting in which the joint effects of spacing, sample size and bandwidth on segmentation can be isolated.

Table 1.

Piecewise-constant intensity function for Level A simulations. For each pulse, the table lists the calendar-time start and end (cal BCE) and the event intensity ${\lambda _k}$ . A low background intensity ${\lambda _{{\rm{bg}}}} = 0.1$ applies across the intervals outside the pulses. Level B simulations use analogous 400-yr pulses in calendar years BP.

Table 1. Long description

A table with three rows and four columns. The columns are labeled Pulse ID, t_start (cal BCE), t_end (cal BCE), and λ_k. Row 1: Pulse ID 1, t_start -3700, t_end -3300, λ_k 1.0. Row 2: Pulse ID 2, t_start -2900, t_end -2500, λ_k 1.5. Row 3: Pulse ID 3, t_start -2100, t_end -1700, λ_k 2.0. The table lists the calendar-time start and end for each pulse and the event intensity λ_k. A low background intensity λ_bg = 0.1 applies across the intervals outside the pulses.

Throughout Level A, pulse width is fixed at 400 yr and separation ${\rm{\Delta }}$ is defined as the centre-to-centre distance between adjacent pulses. This distinction matters. When ${\rm{\Delta }} = 400$ yr, adjacent pulse edges abut; when $\Delta {\rm{ }} \lt 400$ yr, adjacent pulses overlap in calendar time; and when $\Delta {\rm{ }} \gt 400$ yr, a trough or lull exists between them. Reviewer concern about the recoverability of 400-year pulses at similar separations is therefore addressed here by making the relation between width and spacing explicit in both the design and the figure captions.

Events are generated by sampling segment membership with probability proportional to intensity $ \times $ duration and drawing times uniformly within segments. This yields true event times $\{ t_i^{{\rm{true}}}\} _{i = 1}^N$ . To mimic posterior uncertainty in calendar time without calibration, each event is perturbed by Gaussian error with $\sigma = 50$ yr, and the perturbed times are treated as posterior draws when constructing KDEs. Figure 1 illustrates the resulting KDEs and persistent minima for a well-separated case and a more challenging configuration.

KDEs are evaluated on a fixed grid of 400 equally spaced calendar points between $ - 4000$ and $ - 1200$ cal BCE. The main design varies the total sample size $N \in \left\{ {50,100,200,500} \right\}$ , the centre-to-centre separation between pulses ${\rm{\Delta }} \in \left\{ {200,400,600,800} \right\}$ yr, and the bandwidth $h \in \left\{ {80,100,120} \right\}$ yr. For each combination of $\left( {N,{\rm{\Delta }},h} \right)$ , 200 datasets are simulated and analysed with the ensemble KDE–minima rule and with a naive single-KDE minima rule that uses minima of a single KDE constructed from a single draw of calendar ages. The bandwidth values are chosen to interrogate the same multi-century regime across all simulations: they are not intended as optimal selectors, but as a plausible decision range for asking whether 400-year pulses separated by several centuries remain recoverable after smoothing. In other words, the goal is conservative segmentation at the scale of broad episodes, not local feature recovery at the finest temporal resolution. Values substantially larger than 120 yr would increasingly force adjacent episodes to merge, whereas substantially smaller values would make spurious local minima more common. The chosen range therefore brackets a practical compromise between resolution and stability while keeping the design tractable (Bronk Ramsey Reference Bronk Ramsey2017; Crema and Bevan Reference Crema and Bevan2021; Silverman Reference Silverman2018; Scott Reference Scott2015; Sheather Reference Sheather2004).

2.3 Level B simulations: calibrated dates and alternative segmentation methods

Level B simulations extend Level A by including radiocarbon measurement error and calibration, and by comparing the ensemble rule with two alternatives: the naive single-KDE minima rule and a Poisson changepoint model fitted to binned counts.

The simulation window spans 2000–5000 cal BP. Three pulses of width 400 yr are centred at 4500, 3500 and 2500 cal BP, with equal intensities and a background intensity ${\lambda _{{\rm{bg}}}} = 0.1$ , so that pulses contribute approximately 87% of expected events. As in Level A, the low background component prevents the simulations from collapsing into an unrealistically clean pulse-only system and approximates the low-rate dated events expected outside major activity clusters. All pulses have the same height in Level B to focus on calibration effects; the structure is analogous to Table 1 but expressed in calendar years BP. As with Level A, this is a deliberately simplified configuration: real sequences will show unequal pulse heights, variable widths and superimposed long-term trends, so the performance estimates reported here should be read as benchmarks under idealised conditions rather than as universal expectations.

For each true event time $t_i^{{\rm{true}}}$ , an expected radiocarbon age ${\mu _i}$ is first obtained from a specified calibration curve by uncalibration. Three calibration setups are considered:

1. 1. Pure IntCal20 (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards and Friedrich2020),

2. 2. Pure SHCal20 (Hogg et al. Reference Hogg, Heaton, Hua, Palmer, Turney, Southon, Bayliss, Blackwell, Boswijk and Bronk Ramsey2020),

3. 3. A simple “mixed” configuration in which each event is randomly assigned to IntCal20 or SHCal20 with equal probability.

These three setups are included not because large curve-specific differences are expected a priori, but to test whether the broad conclusions of the method are robust to reasonable hemispheric calibration structure. The mixed configuration is therefore a simple heterogeneity stress test rather than a literal demographic scenario.

In all cases, a laboratory measurement ${A_i}$ is then simulated from a normal distribution with mean ${\mu _i}$ and standard deviation ${\sigma _A} = 30$ ${{\rm{\;}}^{14}}$ C yr. These simulated radiocarbon ages and their associated errors are recalibrated using the same curve assignment that generated ${\mu _i}$ , yielding posterior discrete distributions for each date. The KDE ensemble is then constructed as in Section 2.1, using $R = 300$ posterior draws and a fixed grid of 400 calendar points between 2000 and 5000 cal BP. The simulation pipeline is therefore: true calendar ages $ \to $ expected radiocarbon ages under a chosen curve $ \to $ laboratory error $ \to $ recalibration $ \to $ Monte Carlo KDE ensemble.

The Level B design varies $N \in \left\{ {50,100,200,500} \right\}$ , $h \in \left\{ {80,100,120} \right\}$ yr and calibration setup, and performs 200 replicates per combination, per segmentation method (ensemble rule, naive KDE, Poisson changepoint). The Poisson model operates on counts in 50-yr bins and assumes a piecewise-constant intensity; the number and locations of changepoints are estimated with the PELT multiple-changepoint algorithm in the changepoint package (cpt.meanvar, test.stat = “Poisson”) using a Schwarz information criterion (SIC, i.e. BIC) penalty on the log-likelihood (Heaton et al. Reference Heaton, Al-assam and Bard2025; de Navascués et al. Reference de Navascués, Burgarella and Jakobsson2025).

2.4 Performance metrics

For both simulation levels, performance is evaluated by comparing the inferred number of pulses $\hat K$ and their boundaries with the true configuration:

• Pulse-count accuracy. For each replicate, $\hat K$ is compared to the true ${K_{{\rm{true}}}} = 3$ . Results are summarised by the probabilities of correct recovery ( ${p_{{\rm{correct}}}}$ ), under-segmentation ( ${p_{{\rm{under}}}}$ , $\hat K{\rm{ }} \lt {K_{{\rm{true}}}}$ ) and over-segmentation ( ${p_{{\rm{over}}}}$ , $\hat K{\rm{ }} \gt {K_{{\rm{true}}}}$ ).

• Boundary error. Boundary accuracy is summarised in two ways. First, for all replicates, an unconditional mean absolute boundary error (yr) is computed using nearest-boundary matching. This metric remains useful descriptively but becomes difficult to compare across methods when they differ strongly in their tendency to over-segment. Second, for replicates with $\hat K = {K_{{\rm{true}}}}$ , a conditional mean absolute boundary error is computed by ordered one-to-one pairing of sorted estimated and true boundaries. This conditional metric provides the cleaner comparison of temporal localisation once the correct pulse count has been recovered. For replicates with $\hat K \ne {K_{{\rm{true}}}}$ , supplementary tables also report the unconditional error for incorrect- $K$ cases.

These metrics separate recovery of the overall pulse structure from temporal localisation once that structure has been obtained. This dual perspective mirrors recent simulation-based evaluations of radiocarbon-frequency methods and avoids treating boundary precision as commensurate across methods with different segmentation tendencies (Contreras and Meadows Reference Contreras and Meadows2014; Bevan and Crema Reference Bevan and Crema2021; Crema Reference Crema2022). A compact synthesis across information regimes is provided in Table 2, while detailed simulation values appear in Supplementary Tables S1–S3.

Table 2.

Condensed qualitative summary of envelope KDE–minima performance across simulation regimes. Statements are approximate syntheses derived from Supplementary Tables S1–S3 for bandwidths 80–120 yr.

Table 2. Long description

A table comparing envelope rule performance across different scenarios and conditions. The table has three columns: Scenario, Conditions, and Envelope rule performance. It contains four rows with different scenarios: Low information, Medium information (calendar time), Medium information (calibrated), and High information. Each row details the conditions and the corresponding envelope rule performance. For example, under Low information with N less than or equal to 50 and any separation, the performance is frequent under- and over-segmentation with mean boundary errors often exceeding 300 years. Under Medium information (calendar time) with N between 100 and 200 and centre-to-centre separations approaching 800 years, recovery improves in the widest-separation settings with p_correct commonly rising into the 0.5-0.8 range and mean boundary errors around 50-135 years. Under Medium information (calibrated) with N between 100 and 200 and three pulses plus background, calibration reduces recovery and broadens boundaries with p_correct typically in the 0.25-0.55 range and mean boundary errors around 200-300 years. Under High information with N 500, the envelope rule becomes more stable with p_correct commonly in the 0.4-0.7 range and mean boundary errors around 130-240 years.