timeOptBMCMCplot: TimeOptBMCMCplot: Generate summary figures and statistics for...

In astrochron: A Computational Tool for Astrochronology

TimeOptBMCMCplot: Generate summary figures and statistics for TimeOptBMCMC results

Description

TTimeOptBMCMCplot: Generate summary figures and statistics for TimeOptBMCMC results.

Usage

timeOptBMCMCplot(dat,mcmcres,pdfpara,nburnin=NULL,fc_dat=NULL)

Arguments

dat	Stratigraphic series for astrochronologic assessment. First column should be depth or height (in meters), second column should be data value.
mcmcres	Results from timeOptBMCMC
pdfpara	Parameter file from timeOptBMCMC
nburnin	Number or MCMC samples to discard as burn-in interval. If set to NULL, a burn-in detection algorithm is implemented, using the median value from second half of the MCMC chain as the threshold value.
fc_dat	Fourier Coefficients for dat (optional)

Details

The TimeOptBMCMC algorithm (Malinverno & Meyers, 2024) is a Bayesian version of TimeOpt (Meyers, 2015), which evaluates stratigraphic data sets to assess (1) the concentration of spectral power at specified target astronomical periods (spectral fit), and (2) eccentricity amplitude modulations within the precession band (envelope fit). The envelope fit can optionally be omitted ('env' set to FALSE).

For a given cyclostratigraphic data set, TimeOptBMCMC calculates the posterior probability density function (PDF) of the axial precession frequency (k), sedimentation rate (u), Solar System frequencies g_i, and Solar System frequencies s_i, using an adaptive Markov chain Monte Carlo approach. The Bayesian priors for Solar System frequencies g_i and s_i are age-specific, following Table 2 in Hoang et al. (2021). The Bayesian prior for the precession frequency is age-specific, based on a polynomial fit to the model results of Farhat et al. (2022), with uncertainties from Waltham (2015). The Bayesian prior for sedimentation rate is defined as a uniform distribution (between 'sedmin' and 'sedmax').

Function timeOptBMCMCplot generates summary figures and stistics for a given analysis, and is automatically called from function timeOptBMCMC. It can also be called independently; when doing so, be sure to use the prepared data that is returned from timeOptBMCMC (see examples below).

For additional guidance on the application of TimeOptBMCMC, please see Malinverno & Meyers (2024).

References

Farhat, M., Auclair-Desrotour, P., Boue, G., Laskar, J., 2022, The resonant tidal evolution of the Earth-Moon distance: Astronomy & Astrophysics, 665, L1.

Meyers, S.R., 2015, The evaluation of eccentricity-related amplitude modulation and bundling in paleoclimate data: An inverse approach for astrochronologic testing and time scale optimization: Paleoceanography, 30, doi:10.1002/2015PA002850.

Meyers, S.R. and Malinverno, A, 2018, Proterozoic Milankovitch cycles and the history of the solar system: Proceedings of the National Academy of Sciences, www.pnas.org/cgi/doi/10.1073/pnas.1717689115.

Malinverno, A. and Meyers, S.R., 2024, Bayesian estimation of past astronomical frequencies, lunar distance, and length of day from sediment cycles: Geochemistry, Geophysics, Geosystems, 25, e2023GC011176.

Hoang, N.H., Mogavero, F., Laskar, J., 2021, Chaotic diffusion of the fundamental frequencies in the Solar System: Astronomy & Astrophysics, 654, A156.

Waltham, 2015, Milankovitch period uncertainties and their impact on cyclostratigraphy: Journal of Sedimentary Research, 85, 990-998.

Examples

 ## Not run: 
# EXAMPLE 1: Xiamaling Cu/Al
# Obtain the Xiamaling Cu/Al dataset from the Astrochron server
CuAl=getData("Xiamaling-CuAl")

# Isolate interval of interest and interpolate the data to the median sampling interval of 0.012 m.
CuAl_0.012=linterp(iso(CuAl,xmin=263.42,xmax=265.496,genplot=FALSE))

# Run timeOptBMCMC and plot results
res = timeOptBMCMC(CuAl_0.012, sedmin=0.3, sedmax=0.4, age=1400, nsamples=50000)

# If you want to re-plot the results from timeOptBMCMC, note that the Cu/Al data
# was detrended (if selected) and standardized in timeOptBMCMC. Use res$dat, returned
# from timeOptBMCMC
timeOptBMCMCplot(res$dat,res$mcmcres,res$pdfpara)

# EXAMPLE 2: Walvis Ridge a* 
# Obtain the Walvis Ridge a* dataset from the Astrochron server
a=getData("1262-a*")

# Isolate interval of interest and interpolate the data to the median sampling interval of 0.02 m.
a_0.02=linterp(iso(a,xmin=117.36,xmax=139.75))

# Run timeOptBMCMC and plot results
res=timeOptBMCMC(a_0.02, sedmin=1.2, sedmax=1.4, age=55, nsamples=50000)

# Re-plot and return posterior distributions of astronomical periods (including grand cycles)
res2=timeOptBMCMCplot(res$dat,res$mcmcres,res$pdfpara)

# What fraction of the g4-g3 posterior samples exceed the modern grand cycle value of 2.4 Myr?
sum(res2[,16]>2400)/length(res2[,16])
 
## End(Not run)

astrochron documentation built on June 26, 2026, 9:09 a.m.