# Eruption/deposition age estimation demonstration¶

This Jupyter notebook demonstrates the eruption (/deposition) age estimation algorithm of Keller, Schoene, and Samperton (2018) as implemented in Chron.jl.

If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!

Hint: shift-enter to run a single cell, or from the Cell menu select Run All to run the whole file. Any code from this notebook can be copied and pasted into the Julia REPL or a .jl script.

## (1) Load required Julia packages¶

In [20]:
# Load (and install if necessary) the Chron.jl package
try
using Chron
catch
using Pkg
using Chron
end

using Statistics, StatsBase, DelimitedFiles, SpecialFunctions
using KernelDensity: kde
using Plots; gr(); default(fmt = :png)


## (2a) Test Bayesian eruption age estimation with a synthetic dataset¶

#### Generate synthetic zircon dataset, drawing from MeltsVolcanicZirconDistribution¶

In [23]:
dt_sigma = 30; # Timescale relative to analytical uncertainty
N = 42; # Number of zircons

# Draw set of pseudorandom ages from MELTS volcanic zircon distribution,
# with true minimum == 0 and analytical sigma == 1
ages = draw_from_distribution(MeltsVolcanicZirconDistribution,N).*dt_sigma + randn(N);
uncert = ones(N)

# Calculate the weighted mean age
(wx, wsigma, mswd) = awmean(ages, uncert)

h = plot(ylabel="Time before eruption (sigma)", xlabel="N", fg_color_legend=:white, legend=:topleft)
plot!(h,1:length(ages),sort(ages),yerror=uncert*2,seriestype=:scatter, label="Synthetic zircon ages")
plot!(h,collect(xlims()), [0,0], label="Eruption age")
plot!(h,collect(xlims()), [wx,wx], label="Weighted mean, MSWD $(round(mswd,sigdigits=2))")  Out[23]: #### Calculate bootstrapped$\ \mathcal{\vec{f}}_{xtal}(t_r)$¶ In [3]: # Bootstrap the crystallization distribution, # accounting for any expected analytical "tail" beyond eruption/deposition dist = BootstrapCrystDistributionKDE(ages, uncert) dist ./= mean(dist) # Normalize # Plot bootstrapped distribution plot(range(0,1.3,length=length(dist)),dist, label="bootstrapped", ylabel="Probability Density", xlabel="Time before eruption (scaled)", legend=:bottomleft, fg_color_legend=:white) plot!(range(0,1,length=100),MeltsVolcanicZirconDistribution,label="original")  Out[3]: #### Run MCMC to estimate eruption/deposition age distribution of synthetic dataset¶ In [4]: # Configure model nsteps = 200000; # Length of Markov chain burnin = 100000; # Number of steps to discard at beginning of Markov chain # Run MCMC tminDist = metropolis_min(nsteps,dist,ages,uncert); # Print results AgeEst = mean(tminDist[burnin:end]); AgeEst_sigma = std(tminDist[burnin:end]); print("\nEstimated eruption age of synthetic dataset:\n$AgeEst +/- $(2*AgeEst_sigma) Ma (2σ)\n (True synthetic age 0 Ma)") # Plot results h = histogram(tminDist[burnin:end],nbins=50,label="Posterior distribution",xlabel="Eruption Age (Ma)",ylabel="N") plot!(h,[0,0],collect(ylims()),line=(3),label="True (synthetic) age",fg_color_legend=:white) plot!(h,[wx,wx],collect(ylims()),line=(3),label="Weighted mean age",legend=:topright) display(h) sleep(0.5) # (just to make sure this section is finished running)  Estimated eruption age of synthetic dataset: 3.874848093407181 +/- 2.4973592126513786 Ma (2σ) (True synthetic age 0 Ma) InterruptException: Stacktrace: [1] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{Array{Int64,1},Array{Float64,1}}) at /Users/cbkeller/.julia/packages/Plots/5srrj/src/plot.jl:167 [2] plot!(::Plots.Plot{Plots.GRBackend}, ::Array{Int64,1}, ::Vararg{Any,N} where N; kw::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:line, :label, :fg_color_legend),Tuple{Int64,String,Symbol}}}) at /Users/cbkeller/.julia/packages/Plots/5srrj/src/plot.jl:158 [3] top-level scope at In[4]:15 ## (2b) Estimate eruption age for real zircon data¶ The example dataset here is from Wotzlaw et al., 2013 FCT+MLX #### Input dataset (Try pasting in your own data here!)¶ In [5]: # Age and one-sigma uncertainty. ages = [28.196, 28.206, 28.215, 28.224, 28.232, 28.241, 28.246, 28.289, 28.308, 28.332, 28.341, 28.359, 28.379, 28.383, 28.395, 28.4, 28.405, 28.413, 28.415, 28.418, 28.42, 28.422, 28.428, 28.452, 28.454, 28.454, 28.458, 28.468, 28.471, 28.475, 28.482, 28.485, 28.502, 28.52, 28.551, 28.561, 28.565, 28.582, 28.584, 28.586, 28.611, 28.638, 28.655] uncert = [0.019, 0.0155, 0.019, 0.0215, 0.018, 0.023, 0.013, 0.029, 0.0175, 0.0315, 0.0095, 0.0245, 0.0255, 0.0175, 0.0235, 0.014, 0.021, 0.022, 0.0125, 0.0135, 0.016, 0.0195, 0.0175, 0.0125, 0.01, 0.014, 0.015, 0.0205, 0.0155, 0.011, 0.0115, 0.0185, 0.0255, 0.014, 0.0125, 0.013, 0.015, 0.014, 0.012, 0.016, 0.0215, 0.0125, 0.0215] # Sort by age (just to make rank-order plots prettier) t = sortperm(ages) ages = ages[t]; uncert = uncert[t];  #### Calculate bootstrapped$\ \mathcal{\vec{f}}_{xtal}(t_r)$¶ In [6]: # Bootstrap the crystallization distribution, # accounting for any expected analytical "tail" beyond eruption/deposition dist = BootstrapCrystDistributionKDE(ages, uncert) dist ./= mean(dist) # Normalize # Plot bootstrapped distribution plot(range(0,1,length=length(dist)),dist, label="bootstrapped f_xtal", ylabel="Probability Density", xlabel="Time before eruption (unitless)", fg_color_legend=:white)  Out[6]: #### Run MCMC to estimate eruption age¶ In [7]: # Configure model nsteps = 400000; # Length of Markov chain burnin = 150000; # Number of steps to discard at beginning of Markov chain # Run MCMC tminDist = metropolis_min(nsteps,dist,ages,uncert); # Consider also: # tminDist = metropolis_min(nsteps,UniformDistribution,ages,uncert); # tminDist = metropolis_min(nsteps,TriangularDistribution,ages,uncert); # tminDist = metropolis_min(nsteps,HalfNormalDistribution,ages,uncert); # tminDist = metropolis_min(nsteps,MeltsVolcanicZirconDistribution,ages,uncert); # Print results AgeEst = mean(tminDist[burnin:end]); AgeEst_sigma = std(tminDist[burnin:end]); print("\nEstimated eruption age:\n$AgeEst +/- \$(2*AgeEst_sigma) Ma (2σ)\n")

# Plot results
h = histogram(tminDist[burnin:end],nbins=100,label="Posterior distribution",xlabel="Eruption Age (Ma)",ylabel="N",legend=:topleft,fg_color_legend=:white)
# plot!(h,[wx,wx],collect(ylims()),line=(3),label="Weighted mean age",legend=:topleft)
display(h)
sleep(0.5)

Estimated eruption age:
28.183871413905244 +/- 0.03824145453445317 Ma (2σ)

In [8]:
# Plot eruption age estimate relative to rank-order plot of raw zircon ages
h = plot(ylabel="Age (Ma)", xlabel="N", legend=:topleft, fg_color_legend=:white)
plot!(h,1:length(ages),ages,yerror=uncert*2,seriestype=:scatter, markerstrokecolor=:auto, label="Observed ages")
plot!(h,[length(ages)],[AgeEst],yerror=2*AgeEst_sigma, markerstrokecolor=:auto, label="Bayesian eruption age estimate",color=:red)
plot!(h,collect(xlims()),[AgeEst,AgeEst],color=:red, label="")

Out[8]:

In [ ]: