These functions run many simulations on randomly-sampled activity values constrained by the measured activity and estimated background. Excess is calculated by pb210_excess() for each simulation. Prediction results are presented as the median result and are constrained by min (5th percentile) and max (95th percentile) values (instead of quadrature-propagated error like pb210_cic() and pb210_crs()). Note that this may take 10 seconds per 1,000 iterations (depending on hardware).

pb210_cic_monte_carlo(
  cumulative_dry_mass,
  activity,
  background = 0,
  model_top = ~pb210_fit_exponential(..1, ..2),
  decay_constant = pb210_decay_constant(),
  n = 1000,
  sample_activity = pb210_sample_norm,
  sample_background = pb210_sample_norm,
  sample_decay_constant = pb210_sample_norm
)

pb210_crs_monte_carlo(
  cumulative_dry_mass,
  activity,
  background = 0,
  inventory = pb210_inventory_calculator(),
  core_area = pb210_core_area(),
  decay_constant = pb210_decay_constant(),
  n = 1000,
  sample_activity = pb210_sample_norm,
  sample_background = pb210_sample_norm,
  sample_decay_constant = pb210_sample_norm
)

# S3 method for pb210_fit_cic_monte_carlo
predict(object, cumulative_dry_mass = NULL, ...)

# S3 method for pb210_fit_crs_monte_carlo
predict(object, cumulative_dry_mass = NULL, ...)

Arguments

cumulative_dry_mass

The cumulative dry mass of the core (in kg), starting at the surface sample and including all samples in the core. These must be greater than 0 and in increasing order.

activity

A vector of measured lead-210 specific activities (in Bq/kg) and associated error. These can have errors::errors().

background

A vector of estimated background lead-210 specific activity (in Bq/kg) and associated error. These can have errors::errors().

model_top

A fit object, such as one generated by pb210_fit_exponential() or a constant specifying the surface excess. The choice of this value has considerable impact on young dates.

decay_constant

The decay contstant for lead-210 (in 1/years). This is an argument rather than a constant because we have found that different spreadsheets in the wild use different decay constants. See pb210_decay_constant().

n

The number of permutations. The default is 1,000, as Sanchez-Cabeza et al. (2014) found that this was the minimum number of iterations needed for Monte-Carlo uncertainty to converge on the quadrature-propagated uncertainty. In general, Sanchez-Cabeza et al. (2014) used n values from 1,000 to 4,000.

sample_activity, sample_background, sample_decay_constant

Random sampler functions such as pb210_sample_norm() that are called n times for the appropriate argument.

inventory

The cumulative excess lead-210 activity (in Bq), starting at the bottom of the core. By default, this is estimated by the default pb210_inventory_calculator(). If specifying a vector of values, ensure that the surface (0 cumulative mass) value is specified.

core_area

The internal area of the corer (in m^2^). This can be calculated from an internal diameter using pb210_core_area().

object

A fit object generated by pb210_crs() or pb210_cic().

...

Unused.

References

Binford, M.W. 1990. Calculation and uncertainty analysis of ^210^Pb dates for PIRLA project lake sediment cores. Journal of Paleolimnology, 3: 253–267. https://doi.org/10.1007/BF00219461

Sanchez-Cabeza, J.-A., Ruiz-Fernández, A.C., Ontiveros-Cuadras, J.F., Pérez Bernal, L.H., and Olid, C. 2014. Monte Carlo uncertainty calculation of ^210^Pb chronologies and accumulation rates of sediments and peat bogs. Quaternary Geochronology, 23: 80–93. https://doi.org/10.1016/j.quageo.2014.06.002

Examples

# simulate a core core <- pb210_simulate_core() %>% pb210_simulate_counting() # calculate ages using the CRS model crs <- pb210_crs_monte_carlo( pb210_cumulative_mass(core$slice_mass), set_errors( core$activity_estimate, core$activity_se ), n = 100 ) predict(crs)
#> # A tibble: 60 x 24 #> age age_min age_max age_values mar mar_min mar_max mar_values inventory #> <dbl> <dbl> <dbl> <list> <dbl> <dbl> <dbl> <list> <dbl> #> 1 2.38 2.34 2.41 <dbl [100… 0.152 0.150 0.155 <dbl [100… 2991. #> 2 7.36 7.25 7.46 <dbl [100… 0.153 0.150 0.155 <dbl [100… 2561. #> 3 12.4 12.2 12.6 <dbl [100… 0.151 0.149 0.154 <dbl [100… 2189. #> 4 17.5 17.3 17.7 <dbl [100… 0.150 0.148 0.152 <dbl [100… 1866. #> 5 22.7 22.4 22.9 <dbl [100… 0.151 0.148 0.154 <dbl [100… 1589. #> 6 27.8 27.5 28.2 <dbl [100… 0.151 0.147 0.153 <dbl [100… 1352. #> 7 33.1 32.6 33.5 <dbl [100… 0.152 0.148 0.154 <dbl [100… 1149. #> 8 38.3 37.8 38.8 <dbl [100… 0.153 0.150 0.157 <dbl [100… 976. #> 9 43.6 43.0 44.2 <dbl [100… 0.149 0.146 0.153 <dbl [100… 827. #> 10 49.0 48.3 49.7 <dbl [100… 0.149 0.145 0.153 <dbl [100… 698. #> # … with 50 more rows, and 15 more variables: inventory_min <dbl>, #> # inventory_max <dbl>, inventory_values <list>, excess <dbl>, #> # excess_min <dbl>, excess_max <dbl>, excess_values <list>, activity <dbl>, #> # activity_min <dbl>, activity_max <dbl>, activity_values <list>, #> # background <dbl>, background_min <dbl>, background_max <dbl>, #> # background_values <list>