gibbs_ad: Gibbs Sampler for Archaeological Dates

In eratosthenes: Archaeological Synchronism

Gibbs Sampler for Archaeological Dates

Description

A Gibbs sampler for dating archaeological events, to fit relative sequences to absolute, calendrical dates, along with rule-based production dates of artifact types. Relative events can be associated with termini post quos (t.p.q.) and termini ante quos (t.a.q.), which are entered as samples from a given probability density function f(t). This function may take any form, a single date (i.e., with a probability of 1), a continuous uniform distribution (any time between two dates), or a bespoke density (as with calibrated radiocarbon dates). Relative events are modeled on a continuous uniform density between the latest antecedent event and earliest subsequent event.

Usage

gibbs_ad(
  sequences,
  finds = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  trim = TRUE,
  rule = "naive"
)

## S3 method for class 'list'
gibbs_ad(
  sequences,
  finds = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  trim = TRUE,
  rule = "naive"
)

Arguments

sequences	A list of relative sequences of elements (e.g., contexts).
finds	Optional. A list of finds related to (contained in) the elements of sequences.
max_samples	Maximum number of samples to run. Default is 10^5.
size	The number of samples to take on each iteration of the main Gibbs sampler. Default is 10^3.
mcse_crit	Criterion for the Monte Carlo standard error to stop the Gibbs sampler, as based on depositional dates and absolute constraints. The number of Monte Carlo samples for production dates is identical to that depositional dates.
tpq	A list containing termini post quos. Each object in the list consists of: id A character ID of the t.p.q., such as a reference or number. assoc The element in code to which the t.p.q. is associated. samples A vector of samples drawn from the appertaining probability density function of that t.p.q.
taq	A list containing termini ante quos. Each object in the list consists of: id A character ID of the t.a.q., such as a reference or number. assoc The element in code to which the t.p.q. is associated. samples A vector of samples drawn from the appertaining probability density function of that t.p.q.
alpha_	An initial t.p.q. to limit any elements which may occur before the first provided t.p.q. Default is -5000.
omega_	A final t.a.q. to limit any elements which may occur after the after the last provided t.a.q. Default is 1950.
trim	A logical value to determine whether elements that occur before the first t.p.q. and after the last t.a.q. should be omitted from the results (i.e., to "trim" elements at the ends of the sequence, whose marginal densities depend on the selection of alpha_ and omega_). Default is TRUE.
rule	The rule for computing an estimated date of production of a find-type, either "earliest", selecting a production date between the earliest deposition of that type and the next most earliest context, or "naive" (the default), which will select a production date any time between the distribution of that "earliest" date and the depositional date of that artifact.

Details

Gibbs sampling is a conventional method for calibrating and estimating radiocarbon dates in light of absolute constraints and relative sequences: see \insertCitebuck_bayesian_1996,buck_bcal_1999,bronk_ramsey_bayesian_2009;textualeratosthenes, the latter of which uses a mixture of Metropolis-Hastings and Gibbs.

In this implementation, two phases of Gibbs sampling are performed: an initial phase for selecting starting values and then the main sampler, with convergence evaluated using Monte Carlo standard errors (MCSE).

The initial Gibbs sampler results in a vector of starting values randomly sampled for each event up to \sqrt{k} runs, where k is the total number of events. Starting values may therefore take some time to assign, but this initial sampling is necessary to avoid a catastrophic collapse due to floating point errors in the initial selection of random values and will also result in closer starting values with respect to marginal densities.

The main Gibbs sampler uses consistent batch means (CBM) determine convergence and hence when to end the main sampling run: there is no motivation to remove burn-in from the main sampling run nor to run multiple chains. CBM is assured to converge in distribution, see \insertCitejones_fixed-width_2006,flegal_markov_2008;textualeratosthenes. A stopping point for the main sampler is therefore determined using the mean of the Monte Carlo standard errors (MCSE) across all random variates, which is the input of mcse_crit (the mean MCSE for all events). The input max_samples indicates the maximum number of simulations to run, but the sampler will stop if the specified criterion of mcse_crit is passed. The default mean MCSE is set at mcse_crit = 0.5, as the MCSE is measured in years (i.e. to allow for an error +/- 1 year), but, to be sure, individual events will have higher or lower MCSE than this mean criterion, whose primary purpose is as a stopping rule.

Note that the MCSE criterion is applied as a stopping rule for depositional dates and external constraints. The number of Monte Carlo samples for production dates of types is chosen to be identical to that need to pass mcse_crit, such that ultimately the final mean MCSE of all variates may differ from that of the criterion. Depending on the conditional structure of the relative sequences and the timescale of investigation, higher or lower MCSE may be more desirable or acceptable.

For the use dates of artifact type production, use, and deposition, see the gibbs_ad_use function.

Value

A list object of class marginals which contains the following:

• deposition A list of samples from the marginal density of each context's depositional date.

• externals A list of samples of the marginal density of each constraint (t.p.q. and t.a.q.]), as conditioned upon the occurrence of other depositional

• production If a finds object has been input, samples of the marginal density of the production date of finds types will be included in the output. If types are attested in trimmed contexts,

• mcse The Monte Carlo standard errors (MCSE) of the random variates (fixed t.p./a.q. will have a MCSE of 0.)

References

\insertAllCited

Examples

x <- c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J")
y <- c("B", "D", "G", "H", "K")
z <- c("F", "K", "L", "M")
contexts <- list(x, y, z)

f1 <- list(id = "find01", assoc = "D", type = c("type1", "form1"))
f2 <- list(id = "find02", assoc = "E", type = c("type1", "form2"))
f3 <- list(id = "find03", assoc = "G", type = c("type1", "form1"))
f4 <- list(id = "find04", assoc = "H", type = c("type2", "form1"))
f5 <- list(id = "find05", assoc = "I", type = "type2")
f6 <- list(id = "find06", assoc = "H", type = NULL)

artifacts <- list(f1, f2, f3, f4, f5, f6)
 
# external constraints
coin1 <- list(id = "coin1", assoc = "B", type = NULL, samples = runif(100,-320,-300))
coin2 <- list(id = "coin2", assoc = "G", type = NULL, samples = seq(37, 41, length = 100))
  # seq(37, 41, length = 100) is equivalent in concept to runif(100, 37, 41)) 
destr <- list(id = "destr", assoc = "J", type = NULL, samples = 79)

tpq_info <- list(coin1, coin2)
taq_info <- list(destr)

result <- gibbs_ad(contexts, finds = artifacts, tpq = tpq_info, taq = taq_info)

eratosthenes documentation built on June 28, 2025, 1:08 a.m.