msd: Mean Squared Displacement of Events

In eratosthenes: Archaeological Synchronism

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.

Usage

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

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

Arguments

marginalized	The results of gibbs_ad.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. See gibbs_ad for details.

Details

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

• Each event is omitted from all relative and absolute sequences, and the function gibbs_ad is re-run to compute a "jackknifed" Monte Carlo mean for that event.

  • The squared difference of this jackknifed Monte Carlo mean and the original is then computed as its squared "displacement" in time.

  • The mean of the squared displacements of all events is then computed and attributed to the omitted event.

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_, 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 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 <- 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))
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)

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

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