Package {eratosthenes}

Title: Archaeological Synchronism
Version: 1.0.2
Description: Estimation of unknown historical or archaeological dates subject to relationships with other relative dates and absolute constraints, derived as marginal densities from the full joint conditional, using a two-stage Gibbs sampler with consistent batch means to assess convergence. Features reporting on Monte Carlo standard errors, as well as tools for rule-based estimation of dates of production and use of artifact types, aligning and checking relative sequences, and evaluating the impact of the omission of relative/absolute events upon one another.
License: GPL (≥ 3)
Imports: stats, graphics, grDevices, Rcpp, Rdpack, paletteer
RdMacros: Rdpack
Encoding: UTF-8
RoxygenNote: 7.3.2
LinkingTo: Rcpp
Suggests: knitr, rmarkdown, testthat (≥ 3.0.0)
VignetteBuilder: knitr
Config/testthat/edition: 3
NeedsCompilation: yes
Packaged: 2026-10-01 21:53:21 UTC; archaeologus
Author: Stephen A. Collins-Elliott ORCID iD [aut, cre]
Maintainer: Stephen A. Collins-Elliott <sce@utk.edu>
Repository: CRAN
Date/Publication: 2026-10-01 22:40:02 UTC

Create an Absolute Constraint Object

Description

Analogous to the list function, to create absolute constraint (terminus post quem or ante quem). The constraint must contain named elements of "id", "assoc", and "samples", with an option to indicate "type".

Usage

absolute(id, assoc, type = NULL, samples)

## S3 method for class 'character'
absolute(id, assoc, type = NULL, samples)

Arguments

id

a character object, giving a unique ID of the constraint

assoc

the element within an events object to which the constraint is associated

type

(optional) a character object, giving the type of constraint (e.g., a ceramic type, coin, radiocarbon dat). Mutiple types/subtypes/classes can be given as a vector. Default is NULL.

samples

a vector of samples drawn from the appertaining probability density function of that constraint

Value

An absolute object.

Examples

# external constraints
coin1 <- absolute(id = "coin1", assoc = "B", samples = runif(100,-320,-300))
coin2 <- absolute(id = "coin2", assoc = "G", type = "RIC2 57", samples = seq(37, 41, length = 100))
  # seq(37, 41, length = 100) is equivalent in concept to runif(100, 37, 41))
destr <- absolute(id = "destr", assoc = "J", samples = 79)

coin1
coin2
destr

Create an Assemblage Object

Description

Analogous to the sequences or constraints function for relative and absolute events, this function collects one or more finds objects into a single object, for input into gibbs_ad_type.

Usage

assemblage(...)

## S3 method for class 'finds'
assemblage(...)

## S3 method for class 'list'
assemblage(...)

Arguments

...

one or more finds objects.

Value

An assemblage object.

Examples

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

finds_all <- assemblage(f1, f2, f3, f4, f5, f6)
finds_all

Create an Constraints Object

Description

Analogous to the sequences function for relative events, this function collects one or more absolute objects into a single object, for input into gibbs_ad.

Usage

constraints(...)

## S3 method for class 'absolute'
constraints(...)

## S3 method for class 'list'
constraints(...)

Arguments

...

one or more absolute objects, or a list of absolute objects.

Value

A constraints object.

Examples

# external constraints
coin1 <- absolute(id = "coin1", assoc = "B", samples = runif(100,-320,-300))
coin2 <- absolute(id = "coin2", assoc = "G", type = "RIC2 57", samples = seq(37, 41, length = 100))
  # seq(37, 41, length = 100) is equivalent in concept to runif(100, 37, 41))
destr <- absolute(id = "destr", assoc = "J", samples = 79)

tpq <- constraints(coin1, coin2)
taq <- constraints(destr)

tpq
taq

Create an Events Object

Description

Analogous to the c function, to create a sequence of unique events as a vector. Elements may not contain names of "alpha" or "omega", which are restricted for quae_antea and quae_postea.

Usage

events(...)

## S3 method for class 'character'
events(...)

Arguments

...

Comma separated character elements, in order from left (earliest) to right (latest).

Value

An events object.

Examples

# "A" before "B", "B" before "C"
x <- events("A", "B", "C")
x

Create an Finds Object

Description

Analogous to the list function, to create an object of a find (e.g., artifact or other element) related to a particular context or event. The find must contain named elements of "id" and "assoc", with optional inputs of "type" and "residual", to be collected into a single object via the assemblage function. Finds which are datable to absolute constraints should be created as an absolute object.

Usage

finds(id, assoc, type = NULL, residual = FALSE)

## S3 method for class 'character'
finds(id, assoc, type = NULL, residual = FALSE)

Arguments

id

a character object, giving a unique ID of the find.

assoc

the element (e.g., context) within an events object to which the find is associated.

type

(optional) a character object, giving the type of constraint (e.g., a ceramic type, coin, radiocarbon dat). Mutiple types/subtypes/classes can be given as a vector. Default is NULL.

residual

(optional) if TRUE, indicates that the object is residual to its associated event (assoc), e.g., had a final deposition to be regarded prior to its context. Supplying residual = TRUE will suppress it from the estimation of production, use, and depositional dates in the function gibbs_ad_type. Default is FALSE.

Value

A finds object.

Examples

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

f1
f2
f3
f4
f5
f6

Gibbs Sampler for Archaeological Dates

Description

A Gibbs sampler for dating archaeological events, to fit relative sequences to absolute, calendrical dates. 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,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  trim = TRUE,
  quiet = FALSE
)

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

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

Arguments

sequences

A sequences object of relative sequences of elements (e.g., contexts).

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 constraints object containing all termini post quos.

taq

A constraints object containing all termini ante quos.

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.

quiet

Whether to supress messages/progress output while function is running. Default is FALSE.

Details

Gibbs sampling is a conventional method for calibrating and estimating radiocarbon dates in light of absolute constraints and relative sequences: see Buck et al. (1996); Buck et al. (1999); Bronk Ramsey (2009).

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 Jones et al. (2006); Flegal et al. (2008). 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_type function.

Value

A list object of class marginals which contains the following:

References

Bronk Ramsey C (2009). “Bayesian Analysis of Radiocarbon Dates.” Radiocarbon, 51, 337–360. doi:10.1017/s0033822200033865.

Buck CE, Cavanagh WG, Litton CD (1996). Bayesian Approach to Interpreting Archaeological Data. John Wiley and Sons, Chichester.

Buck CE, Christen JA, James GN (1999). “BCal: An On-line Bayesian Radiocarbon Calibration Tool.” Internet Archaeology, 7. doi:10.11141/ia.7.1, https://intarch.ac.uk/journal/issue7/buck/.

Flegal JM, Haran M, Jones GL (2008). “Markov Chain Monte Carlo: Can We Trust the Third Significant Figure?” Statistical Science, 23, 250–260. doi:10.1214/08-STS257.

Jones GL, Haran M, Caffo BS, Neath R (2006). “Fixed-Width Output Analysis for Markov Chain Monte Carlo.” Journal of the American Statistical Association, 101, 1537–1547. doi:10.1198/016214506000000492.

Examples

x <- events("A", "B", "C", "D", "E", "F", "G", "H", "I", "J")
y <- events("B", "D", "G", "H", "K")
z <- events("F", "K", "L", "M")
contexts <- sequences(x, y, z)
 
# external constraints
coin1 <- absolute(id = "coin1", assoc = "B", type = NULL, samples = runif(100,-320,-300))
coin2 <- absolute(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 <- absolute(id = "destr", assoc = "J", type = NULL, samples = 79)

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

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

Gibbs Sampler for Archaeological Dates: Artifact Types

Description

Estimate a densities for the production, use, and deposition dates of an artifact or artifact type. Multiple artifacts and types can be given, which will be pooled into a single type. For example, one can input several individual finds via their id number as comprising a type, or multiple (sub)types/classes as a single type, (e.g., "MGS V amphora", "MGS VI amphora", "MGS V/VI amphora" to construct one group). Depending on whether one is using id numbers or type(s), the id or type argument is used, which takes a vector of the entries' names. The gibbs_ad_type function works on the basis of the presence/absence of types in contexts, sampling a use date between the production and deposition. The stipulation of a rule to determine production dates (naive or earliest) is required.

Usage

gibbs_ad_type(
  sequences,
  finds = NULL,
  id = NULL,
  type = NULL,
  type_name = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  trim = TRUE,
  rule = "naive",
  quiet = FALSE
)

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

Arguments

sequences

A sequences object of relative sequences of elements (e.g., contexts).

finds

An assemblage object of finds related to (contained in) the elements of sequences.

id

A vector of the id of one or more specific finds whose use date is to be estimated. The values of id must match those in the list of finds. If type is used, id is ignored.

type

A vector of one or more types to estimate a use density for. Must contain a value if id is NULL.

type_name

A customized label for the type (e.g., if one is selecting via id or has combined subtypes). If only type is used to select finds, the default will be that label Otherwise the default is simply "Type."

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. Only the MCSE of the use date is used as a stopping rule.

tpq

A constraints object containing all termini post quos.

taq

A constraints object containing all termini ante quos.

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.

quiet

Whether to supress messages/progress output while function is running. Default is FALSE.

Details

See gibbs_ad for information on consistent batch means and Monte Carlo standard error, which are used to determined convergence for the use date.

Value

A list of class type_marginals of the density of a use date, conditional upon production and depositional dates.

Examples

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

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

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

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

# use dates by specifying ids
gibbs_ad_type(contexts, artifacts, id = c("find04", "find05"),
              max_samples = 2000, mcse_crit = 2, tpq = tpq_info, taq = taq_info)

# use dates by specifying types
gibbs_ad_type(contexts, artifacts, type = "type1",
              max_samples = 2000, mcse_crit = 2, tpq = tpq_info, taq = taq_info)

Histogram of Marginal Densities

Description

Wrapper around hist to plot density histograms for select marginal densities (up to 12) in a single plot, from the results of gibbs_ad, or to plot density histograms of the production, deposition, and use of a type, from the results of gibbs_ad_type.

Usage

histogram(
  x,
  events = NULL,
  aspect = c("production", "use", "deposition"),
  breaks = "Freedman-Diaconis",
  xlim = NULL,
  ylim = NULL,
  xlab = "Year",
  palette = NULL,
  opacity = 1,
  legend_pos = "topright"
)

## S3 method for class 'marginals'
histogram(
  x,
  events = NULL,
  aspect = NULL,
  breaks = "Freedman-Diaconis",
  xlim = NULL,
  ylim = NULL,
  xlab = "Year",
  palette = NULL,
  opacity = 1,
  legend_pos = "topright"
)

## S3 method for class 'type_marginals'
histogram(
  x,
  events = NULL,
  aspect = c("production", "use", "deposition"),
  breaks = "Freedman-Diaconis",
  xlim = NULL,
  ylim = NULL,
  xlab = "Year",
  palette = NULL,
  opacity = 0.5,
  legend_pos = "topright"
)

Arguments

x

A list object of class marginals, the output of gibbs_ad, or of class type_marginals, to plot the output of gibbs_ad_type].

events

If plotting a marginals object, a vector or element of the event names to plot. Maximum number of events is 12.

aspect

If plotting a type_marginals object, that is, the output of gibbs_ad_type, a vector of one or more of c("production", "use", "deposition"). The default is all three.

breaks

The number or method of breaks in the histogram. Default is "Freedman-Diaconis". See hist for more.

xlim

The limits of the x-axis. Default is set to the min/max values of all samples.

ylim

The limits of the y-axis. This may need to be adjusted if densities have an extremely narrow interval.

xlab

Label for the x-axis. Default is "Year".

palette

A vector providing the color palette of the histogram. The default is "colorBlindness::paletteMartin" (see palettes_d).

opacity

The opacity/transparency of the histograms for visualizing overlapping events, a value between 0 and 1 (default).

legend_pos

The position of the legend in the plot. Default is "topright".

Details

See also also tidy_marginals for exporting the results of these functions into tidy data frame for custom plotting in e.g., ggplot2.

Value

A density histogram of the selected events/aspects.

A density histogram of the selected events/aspect.

Examples

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

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

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

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

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

# deposition of "B"
histogram(result, "B")

# deposition of "coin2" and deposition of "G"
histogram(result, c("coin2", "G"), opacity = 0.5)

# production, use, and deposition of "type1"
result_type1 <- gibbs_ad_type(contexts, artifacts, type = "type1",
                          max_samples = 3000, mcse_crit = 2)
histogram(result_type1)

Ids of Types

Description

Given a list object of finds (with keys of id, assoc, type in each entry), return a vector of the id elements that belong to one or more specified type.

Usage

ids_of_types(input, type = NULL)

## S3 method for class 'assemblage'
ids_of_types(input, type = NULL)

Arguments

input

An assemblage, comprising finds.

type

A character vector or element

Value

A character vector of ids within a list object of finds class,

Examples

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

artifacts <- assemblage(f1, f2, f3, f4, f5, f6)

ids_of_types(artifacts, type = "type1")
ids_of_types(artifacts, type = c("type1", "type2"))

Mean Squared Displacement of Events

Description

Computes the mean squared displacement (MSD) of all events contained in the relative sequences and absolute constraints used in the execution of gibbs_ad. MSD is not intended for finds, in their production, use, and depositional dates, since these aspects are themselves contingent upon the variates of relative/absolute events.

Usage

msd(
  marginalized,
  sequences,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  quiet = FALSE
)

## S3 method for class 'marginals'
msd(
  marginalized,
  sequences,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  quiet = FALSE
)

Arguments

marginalized

An object of class marginals, the output of gibbs_ad.

sequences

A sequences object of relative sequences of elements (e.g., contexts) used to compute marginalized.

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. A higher MCSE is recommended for situations with a higher number of events in order to reduce computational time.

tpq

A list containing termini post quos used to compute marginalized. See gibbs_ad for details.

taq

A list containing termini ante quos used to compute marginalized. See gibbs_ad for details.

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.

quiet

Whether to supress messages/progress output while function is running. Default is FALSE.

Details

The MSD entails the following jackknife/leave-one-out style routine:

If an event has a low MSD, it bears a low impact on the rest of the events within the full joint conditional density. If it is has a high MSD, other events depend heavily upon its inclusion in the full joint density.

Trimming is not implemented in the computation of MSD, and so attention should be paid to the selection of alpha_ and omega_, which should be reported. This is owing to the way in which, if an absolute constraint (tpq or taq) is omitted that happens to be an earliest or latest bounding event, there still needs to be earliest and latest thresholds in place.

This function is fairly computationally intensive and thus a lower value of max_samples and a higher value of mcse_crit may be warranted.

Value

Output is a list containing a data frame MSD_stats giving the mean MC date, the MCSE, the MSD, the variance of the squared displacements (not the standard error), and sample size, as well as a vector bounds of the values of alpha_ and omega_.

Examples

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

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

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

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

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

result_msd <- msd(result, contexts, max_samples = 5000,
                  mcse_crit = 2, tpq = tpq_info, taq = taq_info)

Quae Antea

Description

For a list of multiple partial sequences of events objects, generates a list which, for each element, giving the elements that occur before it ("quae antea"). This is analogous to a recursive trace through all partial sequences from right to left. An element "alpha" is added to all sets to avoid empty vectors. See also quae_postea.

Usage

quae_antea(...)

## S3 method for class 'events'
quae_antea(...)

## S3 method for class 'list'
quae_antea(...)

## S3 method for class 'sequences'
quae_antea(...)

Arguments

...

Objects of class events, or a sequences object, a valid list of events.

Value

A list of vector objects, which contain the elements that occur before any one given element in the input sequences.

Examples

x <- events("A", "B", "C")
y <- events("B", "D", "E", "C", "F")
z <- events("C", "G")

quae_antea(x, y, z)

Quae Postea

Description

For a list of multiple partial sequences (of vector objects), generate another list which, for each element, gives all elements that occur after it ("quae postea"). This is analogous to a recursive trace through all partial sequences from left to right. A final element "omega" is added to all sets to avoid empty vectors. See also quae_antea.

Usage

quae_postea(...)

## S3 method for class 'events'
quae_postea(...)

## S3 method for class 'list'
quae_postea(...)

## S3 method for class 'sequences'
quae_postea(...)

Arguments

...

Objects of class events, or a sequences object, a valid list of events.

Value

A list of vector objects, which contain the elements that occur after any one given element in the input sequences.

Examples

x <- events("A", "B", "C")
y <- events("B", "D", "E", "C", "F")
z <- events("C", "G")

quae_postea(x)
quae_postea(x, y, z)

a <- sequences(x, y, z)
quae_postea(a)

Adjust Sequence to Target

Description

Given an "input" sequence of elements and another "target" sequence that contains fewer elements in a different order, shift the order of the input sequence to match that of the target, keeping all other elements as proximate to one another as possible. This adjusted ranking is accomplished using piecewise linear interpolation between joint elements ranks. That is, joint rankings are plotted, with input rankings along the x axis and target rankings on the y axis. Remaining rankings in the input sequence are assigned a ranking of y based on the piecewise linear function between joint rankings. If the rank order of elements in the target are identical to those in the input, the result is identical to the input. A minimum number of three joint elements in both the input and target are required.

Usage

seq_adj(input, target)

## S3 method for class 'events'
seq_adj(input, target)

Arguments

input

An events object of unique ordered elements.

target

An events object of unique ordered elements. containing at least three of the same elements as input.

Value

An events object of the adjusted sequence.

Examples

x <- events("A", "B", "C", "D", "E", "F", "G", "H", "I", "J") # the input sequence of events
y <- events("D", "A", "J") # the target sequence of events

seq_adj(x, y)

Sequence Diagnostic

Description

If the creation of a sequences object has failed, this function checks events objects for instances of disagreement, by proceeding through all events in order, agglomerating them and checking for sequence validty. Hence, there are events which are assured to be valid, these should be placed first in the input to seq_diag. If there is no information about the validity of the events, the shuffle option can be set to TRUE, which will randomly permute the order in which events are agglomerated.

Usage

seq_diag(..., shuffle = FALSE)

## S3 method for class 'events'
seq_diag(..., shuffle = FALSE)

## S3 method for class 'list'
seq_diag(..., shuffle = FALSE)

Arguments

...

Objects of events class.

shuffle

Whether to randomly permute the order of the events. Default is FALSE

Value

A sequences object.

Examples

u <- events("A", "D", "E")
v <- events("E", "D")
w <- events("B", "F", "C")
x <- events("A", "B", "C", "D", "E")

seq_diag(u, v, w, x)

a <- list(u, v, w, x)
seq_diag(a)

Create a Sequences Object

Description

Analogous to the list function, a sequences object contains multiple events objects.

Usage

sequences(...)

## S3 method for class 'events'
sequences(...)

## S3 method for class 'list'
sequences(...)

Arguments

...

objects of events class, or a list of events objects.

Value

A sequences object.

Examples

x <- events("A", "B", "C", "D", "E")
y <- events("B", "D", "F")
z <- events("A", "C", "F", "G")
sequences(x, y, z)

a <- list(x, y, z)
sequences(a)

Squared Displacement for a Target Event

Description

Computes the squared displacement for a target event within the joint conditional density, estimating how much the omission of every other event will change the date of the target. See also msd. If the target event is a find or type, the displacement of the use date is used, since use is contingent upon both production and deposition.

Usage

sq_disp(
  marginalized,
  target = NULL,
  sequences,
  finds = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  rule = "naive",
  quiet = FALSE
)

## S3 method for class 'marginals'
sq_disp(
  marginalized,
  target = NULL,
  sequences,
  finds = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  rule = NULL,
  quiet = FALSE
)

## S3 method for class 'type_marginals'
sq_disp(
  marginalized,
  target = NULL,
  sequences,
  finds = NULL,
  max_samples = 10^5,
  size = 10^3,
  mcse_crit = 0.5,
  tpq = NULL,
  taq = NULL,
  alpha_ = -5000,
  omega_ = 1950,
  rule = "naive",
  quiet = FALSE
)

Arguments

marginalized

The results of gibbs_ad or gibbs_ad_type.

target

The target event (any event for which to estimate squared displacement. If using the results of gibbs_ad_type, that type is by default the target (otherwise, for sequences/t.p.q./t.a.q. one should use the output of gibbs_ad).

sequences

A list of relative sequences of elements (e.g., contexts) used to compute marginalized.

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. A higher MCSE is recommended for situations with a higher number of events in order to reduce computational time.

tpq

A list containing termini post quos used to compute marginalized. See gibbs_ad for details.

taq

A list containing termini ante quos used to compute marginalized. See gibbs_ad for details.

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.

rule

The rule for computing an estimated date of production, if using an artifact type as a target date. See gibbs_ad_type for details.

quiet

Whether to supress messages/progress output while function is running. Default is FALSE.

Details

Displacement is computed via the following jackknife/leave-one-out-style routine:

If an event has a low squared displacement, it has a low impact on the dating of the target event. If it is has a high squared displacement, the target event's date depends heavily upon its inclusion in the full joint density.

Trimming is not implemented in the estimation of squared displacement, and so attention should be paid to the selection of alpha_ and omega_, and reported. This is owing to the way in which, if an absolute constraint (tpq or taq) is omitted that happens to be an earliest or latest bounding event, there still needs to be earliest and latest thresholds in place.

This function is fairly computationally intensive, and so a lower value of max_samples or higher value of mcse_crit may be warranted.

Value

Output is a list containing a data frame sq_disp giving the displacement with respect to all other events and a vector bounds of the values of alpha_ and omega_.

Examples

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

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

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

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

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

# max_samples lowered and msce_crit raised for examples

# squared displacement for depositional context "E"
E_sqdisp <- sq_disp(result, target = "E", sequences = contexts, 
                    max_samples = 3000, mcse_crit = 2, tpq = tpq_info, taq = taq_info)

result_type1 <- gibbs_ad_type(contexts, finds = artifacts, type = "type1",
                              tpq = tpq_info, taq = taq_info)

# squared displacement for production of artifact type "type1"
type1_sqdisp <- sq_disp(result_type1, sequences = contexts, finds = artifacts,
                        max_samples = 3000, mcse_crit = 2, tpq = tpq_info, taq = taq_info)

Synthetic Ranking

Description

Using a sequences object of two or more partial sequences, all of which observe the same order of elements, create a single "synthetic" ranking. This is accomplished by counting the total number of elements after running a recursive trace through all partial sequences (via quae_postea). If partial sequences are inconsistent in their rankings, a NULL value is returned.

Usage

synth_rank(obj, ties = "average")

## S3 method for class 'sequences'
synth_rank(obj, ties = "average")

Arguments

obj

A sequences object.

ties

The way in which ties are handled per the rank function. The default is "ties = average".

Value

An events object containing the synthesized ranking.

Examples

x <- events("A", "B", "C", "D", "E")
y <- events("B", "D", "F", "E")
a <- sequences(x, y)

synth_rank(a)

Convert Marginals to Tidy (Molten) Data Frame

Description

Takes the results of gibbs_ad or gibbs_ad_type and "melts" the list into a tidy data frame (Wickham 2014). Each row of the molten data frame will contain the index of the Monte Carlo sample, the sample itself, and then the event name.

Usage

tidy_marginals(input)

## S3 method for class 'marginals'
tidy_marginals(input)

## S3 method for class 'type_marginals'
tidy_marginals(input)

Arguments

input

An object of class marginals or type_marginals, the output of gibbs_ad or gibbs_ad_type.

Value

A data frame giving the MC sampling index (idx), the sample (year), and the event (event).

References

Wickham H (2014). “Tidy Data.” Journal of Statistical Software, 59, 1–23. doi:10.18637/jss.v059.i10.

Examples

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

# external constraints
coin1 <- absolute(id = "coin1", assoc = "B", type = NULL, samples = runif(100,-320,-300))
coin2 <- absolute(id = "coin2", assoc = "G", type = NULL, samples = seq(37, 41, length = 100))
destr <- absolute(id = "destr", assoc = "J", type = NULL, samples = 79)

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

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

tidy_marginals(result)

Traceplot of Gibbs Samples

Description

Wrapper around plot to make a traceplot of Gibbs samples from gibbs_ad. See histogram. Maximum number of simultaneous events to display is 12. for plotting a density histogram of events.

Usage

traceplot(
  x,
  events = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Index",
  ylab = "Year",
  palette = NULL,
  opacity = 1,
  legend_pos = "topright"
)

## S3 method for class 'marginals'
traceplot(
  x,
  events = NULL,
  xlim = NULL,
  ylim = NULL,
  xlab = "Index",
  ylab = "Year",
  palette = NULL,
  opacity = 1,
  legend_pos = "topright"
)

Arguments

x

A list object of class marginals, the output of gibbs_ad.

events

A vector or element of the event names to plot. Maximum number of events is 12.

xlim

The limits of the x-axis (optional).

ylim

The limits of the y-axis (optional).

xlab

Label for the y-axis. Default is "Index".

ylab

Label for the y-axis. Default is "Year".

palette

A vector providing the color palette of the histogram. The default is "colorBlindness::paletteMartin" (see palettes_d).

opacity

The opacity/transparency of the traceplot, if visualizing overlapping events. A value between 0 and 1 (default).

legend_pos

The position of the legend in the plot. Default is "topright".

Details

See also tidy_marginals for exporting the results of these functions into tidy data frame for custom plotting in e.g., ggplot2.

Value

A traceplot of the Gibbs samples of the selected events.

Examples

x <- events("A", "B", "C", "D", "E", "F", "G", "H", "I", "J")
y <- events("B", "D", "G", "H", "K")
z <- events("F", "K", "L", "M")
contexts <- sequences(x, y, z)
 
# external constraints
coin1 <- absolute(id = "coin1", assoc = "B", type = NULL, samples = runif(100,-320,-300))
coin2 <- absolute(id = "coin2", assoc = "G", type = NULL, samples = seq(37, 41, length = 100))
destr <- absolute(id = "destr", assoc = "J", type = NULL, samples = 79)

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

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

traceplot(result, "B")
traceplot(result, c("coin1", "B", "H"), opacity = 0.5)