{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bayesian eruption age estimation example notebook\n",
"\n",
"This Jupyter notebook demonstrates the Bayesian eruption age estimation approach of [Keller, Schoene, and Samperton (2018)](https://doi.org/10.7185/geochemlet.1826)\n",
"\n",
"For future work, see [Chron.jl](https://github.com/brenhinkeller/Chron.jl), which can be run in [standalone eruption/deposition mode](https://github.com/brenhinkeller/Chron.jl#standalone-age-depth-modelling) similar to this notebook.\n",
"\n",
"Hint: `shift`-`enter` to run a single cell, or from the menu select `Cell` > `Run All` to run the whole file\n",
"***\n",
"## (1) Load external resources \n",
"#### (run this first!)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Required packages\n",
"using Plots; gr();\n",
"using KernelDensity: kde\n",
"using Statistics, SpecialFunctions\n",
"\n",
"# Download file from the BayeZirChron repo if necessary\n",
"if ~isfile(\"DistTools.jl\")\n",
" download(\"https://raw.githubusercontent.com/brenhinkeller/BayeZirChron.c/master/julia/DistTools.jl\",\"DistTools.jl\")\n",
"end\n",
"\n",
"# Functions we'll be using here\n",
"include(\"DistTools.jl\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## (2a) Test bayesian eruption age estimation with a synthetic dataset\n",
"#### Generate synthetic zircon dataset, drawing from `MeltsVolcanicZirconDistribution`"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dt_sigma = 30 # Timescale relative to analytical uncertainty\n",
"N = 42 # Number of zircons\n",
" \n",
"# Draw set of pseudorandom ages from MELTS volcanic zircon distribution, \n",
"# with true minimum == 0 and analytical sigma == 1\n",
"ages = draw_from_distribution(MeltsVolcanicZirconDistribution,N).*dt_sigma + randn(N)\n",
"uncert = ones(N)\n",
"\n",
"# Calculate the weighted mean age\n",
"(wx, wsigma, mswd) = awmean(ages, uncert)\n",
"\n",
"h = plot(ylabel=\"Time before eruption (sigma)\", xlabel=\"N\", fg_color_legend=:white, legend=:topleft)\n",
"plot!(h,1:length(ages),sort(ages),yerror=uncert*2,seriestype=:scatter, markerstrokecolor=:auto, label=\"Synthetic zircon ages\")\n",
"plot!(h,collect(xlims()), [0,0], label=\"Eruption age\")\n",
"plot!(h,collect(xlims()), [wx,wx], label=\"Weighted mean, MSWD $(round(mswd,digits=2))\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Calculate bootstrapped $\\ \\mathcal{\\vec{f}}_{xtal}(t_r)$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Maximum extent of expected analytical tail (beyond eruption/deposition)\n",
"maxTailLength = mean(uncert) * norm_quantile(1 - 1/(length(ages) + 1));\n",
"included = (ages .- minimum(ages)) .>= maxTailLength;\n",
"\n",
"# Bootstrapped crystallization distribution, excluding maximum analytical tail\n",
"if sum(included) > 5\n",
" dist = BootstrapCrystDistributionKDE(ages[included]);\n",
"else\n",
" # Avoid edge cases at n = 0 and n = 2;\n",
" # Default to n = 1 instead, which yields a half-normal distribution\n",
" dist = BootstrapCrystDistributionKDE([0.]);\n",
"end\n",
"dist = dist./mean(dist);\n",
"\n",
"# Plot bootstrapped distribution\n",
"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)\n",
"plot!(range(0,1,length=100),MeltsVolcanicZirconDistribution,label=\"original\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Run MCMC to estimate eruption/deposition age distribution of synthetic dataset"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Estimated eruption age of synthetic dataset:\n",
" 0.03542759403584769 +/- 2.639939347467269 Ma (2σ)\n",
" (True synthetic age 0 Ma)"
]
}
],
"source": [
"# Configure model\n",
"nsteps = 200000 # Length of Markov chain\n",
"burnin = 100000 # Number of steps to discard at beginning of Markov chain\n",
"\n",
"# Run MCMC\n",
"tminDist = zeros(nsteps)\n",
"metropolis_min(nsteps,dist,ages,uncert,tminDist)\n",
"\n",
"# Print results\n",
"AgeEst = mean(tminDist[burnin:end])\n",
"AgeEst_sigma = std(tminDist[burnin:end])\n",
"print(\"\\nEstimated eruption age of synthetic dataset:\\n $AgeEst +/- $(2*AgeEst_sigma) Ma (2σ)\\n (True synthetic age 0 Ma)\")\n",
"\n",
"# Plot results\n",
"h = histogram(tminDist[burnin:end],nbins=50,label=\"Posterior distribution\",xlabel=\"Eruption Age (Ma)\",ylabel=\"N\")\n",
"plot!(h,[0,0],collect(ylims()),line=(3),label=\"True (synthetic) age\",fg_color_legend=:white)\n",
"plot!(h,[wx,wx],collect(ylims()),line=(3),label=\"Weighted mean age\",legend=:topright)\n",
"display(h)\n",
"\n",
"sleep(0.5) # (just to make sure this section is finished running)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## (2b) Estimate eruption age for real zircon data\n",
"The example dataset here is from [Wotzlaw et al., 2013](https://doi.org/10.1130/G34366.1) FCT+MLX\n",
"#### Input dataset (Try pasting in your own data here!)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Age and one-sigma uncertainty.\n",
"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]\n",
"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]\n",
"\n",
"# Sort by age (just to make rank-order plots prettier)\n",
"t = sortperm(ages)\n",
"ages = ages[t];\n",
"uncert = uncert[t];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Calculate bootstrapped $\\ \\mathcal{\\vec{f}}_{xtal}(t_r)$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Maximum extent of expected analytical tail (beyond eruption/deposition)\n",
"maxTailLength = mean(uncert) * norm_quantile(1-1/(length(ages) + 1))\n",
"included = (ages .- minimum(ages)) .>= maxTailLength\n",
"\n",
"# Bootstrapped crystallization distribution, excluding maximum analytical tail\n",
"if sum(included) > 5\n",
" dist = BootstrapCrystDistributionKDE(ages[included])\n",
"else\n",
" # Avoid edge cases at n = 0 and n = 2;\n",
" # Default to n = 1 instead, which yields a half-normal distribution\n",
" dist = BootstrapCrystDistributionKDE([0.])\n",
"end\n",
"dist = dist./mean(dist)\n",
"\n",
"# Plot bootstrapped distribution\n",
"plot(range(0,1,length=length(dist)),dist, label=\"bootstrapped\", ylabel=\"Probability Density\", xlabel=\"Time before eruption (unitless)\", fg_color_legend=:white)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Run MCMC to estimate eruption age"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Estimated eruption age:\n",
" 28.1839792870061 +/- 0.037925204048165476 Ma (2σ)\n"
]
}
],
"source": [
"# Configure model\n",
"nsteps = 400000 # Length of Markov chain\n",
"burnin = 150000 # Number of steps to discard at beginning of Markov chain\n",
"\n",
"# Run MCMC\n",
"tminDist = zeros(nsteps)\n",
"metropolis_min(nsteps,dist,ages,uncert,tminDist)\n",
"# Consider also:\n",
"# metropolis_min(nsteps,UniformDistribution,ages,uncert,tminDist);\n",
"# metropolis_min(nsteps,TriangularDistribution,ages,uncert,tminDist);\n",
"# metropolis_min(nsteps,HalfNormalDistribution,ages,uncert,tminDist);\n",
"# metropolis_min(nsteps,MeltsVolcanicZirconDistribution,ages,uncert,tminDist);\n",
"\n",
"# Print results\n",
"AgeEst = mean(tminDist[burnin:end]);\n",
"AgeEst_sigma = std(tminDist[burnin:end]);\n",
"print(\"\\nEstimated eruption age:\\n $AgeEst +/- $(2*AgeEst_sigma) Ma (2σ)\\n\")\n",
"\n",
"# Plot results\n",
"h = histogram(tminDist[burnin:end],nbins=100,label=\"Posterior distribution\",xlabel=\"Eruption Age (Ma)\",ylabel=\"N\",legend=:topleft,fg_color_legend=:white)\n",
"# plot!(h,[wx,wx],collect(ylims()),line=(3),label=\"Weighted mean age\",legend=:topleft)\n",
"display(h)\n",
"sleep(0.5)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Plot eruption age estimate relative to rank-order plot of raw zircon ages\n",
"h = plot(ylabel=\"Age (Ma)\", xlabel=\"N\", legend=:topleft, fg_color_legend=:white)\n",
"plot!(h,1:length(ages),ages,yerror=uncert*2,seriestype=:scatter, markerstrokecolor=:auto, label=\"Observed ages\")\n",
"plot!(h,[length(ages)],[AgeEst],yerror=2*AgeEst_sigma, markerstrokecolor=:auto, label=\"Bayesian eruption age estimate\",color=:red)\n",
"plot!(h,collect(xlims()),[AgeEst,AgeEst],color=:red, label=\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.6.3",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}