Isoplot

intcal13
NamedTuple with 4 fields:
  dt            = Int64 1
  age_calendar  = StepRange{Int64, Int64}   1:1:50000
  age_sigma     = Vector{Float64}(50000,)   [8.0 ... 274.0]
  age_14C       = Vector{Float64}(50000,)   [196.8 ... 46401.0]

The IntCal13 radiocarbon calibration curve

Reimer PJ, Bard E, Bayliss A, Beck JW, Blackwell PG, Bronk Ramsey C, Buck CE, Cheng H, Edwards RL, Friedrich M, Grootes PM, Guilderson TP, Haflidason H, , Hajdas I, Hatté C, Heaton TJ, Hogg AG, Hughen KA, Kaiser KF, Kromer B, Manning SW, Niu M, Reimer RW, Richards DA, Scott EM, Southon JR, Turney CSM, van der Plicht J. (2013) "IntCal13 and MARINE13 radiocarbon age calibration curves 0-50000 years calBP". Radiocarbon 55(4). doi: 10.2458/azujsrc.55.16947

intcal20
NamedTuple with 4 fields:
  dt            = Int64 1
  age_calendar  = StepRange{Int64, Int64}   1:1:55000
  age_sigma     = Vector{Float64}(55000,)   [11.0 ... 1024.0]
  age_14C       = Vector{Float64}(55000,)   [197.0 ... 50100.0]

The IntCal20 northern hemisphere radiocarbon calibration curve.

Reimer P, Austin WEN, Bard E, Bayliss A, Blackwell PG, Bronk Ramsey C, Butzin M, Cheng H, Edwards RL, Friedrich M, Grootes PM, Guilderson TP, Hajdas I, Heaton TJ, Hogg AG, Hughen KA, Kromer B, Manning SW, Muscheler R, Palmer JG, Pearson C, van der Plicht J, Reimer RW, Richards DA, Scott EM, Southon JR, Turney CSM, Wacker L, Adolphi F, Büntgen U, Capano M, Fahrni S, Fogtmann-Schulz A, Friedrich R, Köhler P, Kudsk S, Miyake F, Olsen J, Reinig F, Sakamoto M, Sookdeo A, Talamo S. (2020). "The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP)". Radiocarbon 62. doi: 10.1017/RDC.2020.41.

The Marine20 radiocarbon calibration curve.

Heaton TJ, Köhler P, Butzin M, et al. (2020) "Marine20—The Marine Radiocarbon Age Calibration Curve (0–55,000 cal BP)" Radiocarbon 62(4):779-820. doi: 10.1017/RDC.2020.68

shcal20
NamedTuple with 4 fields:
  dt            = Int64 1
  age_calendar  = StepRange{Int64, Int64}   1:1:55000
  age_sigma     = Vector{Float64}(55000,)   [10.0 ... 1008.0]
  age_14C       = Vector{Float64}(55000,)   [171.0 ... 50132.0]

The SHCal20 southern hemisphere radiocarbon calibration curve.

Hogg AG, Heaton TJ, Hua Q, Palmer JG, Turney CSM, Southon J, Bayliss A, Blackwell PG, Boswijk G, Bronk Ramsey C, Pearson C, Petchey F, Reimer P, Reimer R, Wacker L. (2020). "SHCal20 Southern Hemisphere calibration, 0-55,000 years cal BP". Radiocarbon 62. doi: 10.1017/RDC.2020.59

struct UPbAnalysis{T} <: NoInitialDaughterAnalysis{T}

Core type for U-Pb analyses. Wraps an Analysis object which has fields

μ :: SVector{T<:AbstractFloat}
σ :: SVector{T<:AbstractFloat}
Σ :: SMatrix{T<:AbstractFloat}

where μ contains the means

μ = [μ²⁰⁷Pb/²³⁵U, μ²⁰⁶Pb/²³⁸U]

where σ contains the standard deviations

σ = [σ²⁰⁷Pb/²³⁵U, σ²⁰⁶Pb/²³⁸U]

and Σ contains the covariance matrix

Σ = [ σ₇_₅^2  σ₇_₅σ₃_₈
     σ₇_₅σ₃_₈  σ₃_₈^2 ]

If σ is not provided, it will be automatically calculated from Σ, given that σ.^2 = diag(Σ).

UPbAnalysis(r²⁰⁷Pb²³⁵U, σ²⁰⁷Pb²³⁵U, r²⁰⁶Pb²³⁸U, σ²⁰⁶Pb²³⁸U, correlation; T=Float64)

Construct a UPbAnalysis object from individual isotope ratios and (1-sigma!) uncertainties.

Examples

julia> UPbAnalysis(22.6602, 0.0175, 0.40864, 0.00017, 0.83183)
UPbAnalysis{Float64}([22.6602, 0.40864], [0.00030625000000000004 2.4746942500000003e-6; 2.4746942500000003e-6 2.8900000000000004e-8])

struct UThPbAnalysis{T} <: NoInitialDaughterAnalysis{T}

Core type for U-Th-Pb analyses. Wraps an Analysis object which has fields

μ :: SVector{T<:AbstractFloat}
σ :: SVector{T<:AbstractFloat}
Σ :: SMatrix{T<:AbstractFloat}

where μ contains the means

μ = [μ²⁰⁷Pb/²³⁵U, μ²⁰⁶Pb/²³⁸U, μ²⁰⁸Pb/²³²Th]

where σ contains the standard deviations

σ = [σ²⁰⁷Pb²³⁵U, σ²⁰⁶Pb²³⁸U, σ²⁰⁸Pb/²³²Th]

and Σ contains the covariance matrix

Σ = [ σ₇_₅^2  σ₇_₅σ₃_₈ σ₇_₅σ₈_₂
     σ₇_₅σ₃_₈  σ₃_₈^2  σ₃_₈σ₈_₂
     σ₇_₅σ₈_₂ σ₃_₈σ₈_₂  σ₈_₂^2 ]

If σ is not provided, it will be automatically calculated from Σ, given that σ.^2 = diag(Σ).

UThPbAnalysis(μ, Σ)

Construct a UPbAnalysis object from vector of isotope ratios and (1-sigma!) covariance matrix.

calibrate(d::UPbSIMSData, calib::UPbSIMSCalibration, [cyclefilter]; 
    blank64::Number = stacey_kramers(0)[1], 
    blank74::Number = stacey_kramers(0)[2], 
    U58::Number = 1/137.818, 
    baseline::Number = 0,
)

Calibrate one or more U-Pb SIMS analyses d with the calibration calib, optionally filtering by cyclefilter, resulting in a UPbAnalysis for each input UPbSIMSData object.

A Pb-204 based blank subtraction (common-Pb subtraction) is performed, using by default a present-day Stacey-Kramers common Pb composition.

U-235 is estimated based on measured U-238, assuming a 235/238 ratio of U57 – by default 1/137.818 (i.e., Heiss et al. 2012, doi: 10.1126/science.1215507)

A baseline-subtraction (equal for all masses) may optionally be performed given baseline excess counts per minute.

Examples

julia> standards = importsimsdata("./examples/data")
10-element Vector{Isoplot.UPbSIMSData{Float64}}:
 Isoplot.UPbSIMSData{Float64}(...)
 Isoplot.UPbSIMSData{Float64}(...)
 ⋮
 Isoplot.UPbSIMSData{Float64}(...)

julia> standardages = fill(1099., length(standards));

julia> calib = calibration(standards, standardages)
Isoplot.UPbSIMSCalibration{Float64}(...)

julia> calibrate(standards[1], calib)
UPbAnalysis{Float64}(Analysis2D{Float64}([1.8930052834566204, 0.1835092671444941], [0.021732106262365054, 0.0011423135956512402], [0.0004722844425987264 1.4604632531997064e-5; 1.4604632531997064e-5 1.304880350809665e-6]))

calibrate_blockwise(d::UPbSIMSData, calib::UPbSIMSCalibration, [cyclefilter]; 
    blank64::Number = stacey_kramers(0)[1], 
    blank74::Number = stacey_kramers(0)[2], 
    U58::Number = 1/137.818, 
    baseline::Number = 0,
    blocksize::Integer = 3,
)

Calibrate one or more U-Pb SIMS analyses d with the calibration calib, on a block-by-block basis into blocks of length blocksize, optionally filtering by cyclefilter, resulting in vector of UPbAnalysis objects for each UPbSIMSData object.

A Pb-204 based blank subtraction (common-Pb subtraction) is performed, using by default a present-day Stacey-Kramers common Pb composition.

U-235 is estimated based on measured U-238, assuming a 235/238 ratio of U57 – by default 1/137.818 (i.e., Heiss et al. 2012, doi: 10.1126/science.1215507)

A baseline-subtraction (equal for all masses) may optionally be performed given baseline excess counts per minute.

Examples

julia> standards = importsimsdata("./examples/data")
10-element Vector{Isoplot.UPbSIMSData{Float64}}:
 Isoplot.UPbSIMSData{Float64}(...)
 Isoplot.UPbSIMSData{Float64}(...)
 ⋮
 Isoplot.UPbSIMSData{Float64}(...)

julia> standardages = fill(1099., length(standards));

julia> calib = calibration(standards, standardages)
Isoplot.UPbSIMSCalibration{Float64}(...)

julia> calibrate_blockwise(standards[1], calib)
6-element Vector{UPbAnalysis{Float64}}:
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.9522725644621401, 0.19479933826738347], [0.10870961850778198, 0.0033669214494005937], [0.011817781156107494 0.0003491283425207373; 0.0003491283425207373 1.1336160046433794e-5]))
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.970015175457096, 0.18617712628434524], [0.01431116289886847, 0.0023768591910531384], [0.00020480938351794938 -1.1409841282689287e-6; -1.1409841282689287e-6 5.649459614093779e-6]))
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.868382104433831, 0.18258981150055945], [0.03784079639497292, 0.0012781398026132764], [0.0014319258718057954 -4.8258045726099764e-5; -4.8258045726099764e-5 1.633641355024305e-6]))
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.8869190705183156, 0.18091237376063943], [0.038037522681535894, 0.0013370053746898942], [0.0014468531317483576 -1.2603028339239038e-6; -1.2603028339239038e-6 1.7875833719496643e-6]))
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.856977011176083, 0.17882174298755873], [0.07424037726958133, 0.0014834460179062856], [0.005511633617129768 8.764295590103794e-5; 8.764295590103794e-5 2.2006120880420158e-6]))
 UPbAnalysis{Float64}(Analysis2D{Float64}([1.8386867975254575, 0.18063423529792652], [0.05211549498674472, 0.00248005466136501], [0.002716024817713414 9.982927603459965e-5; 9.982927603459965e-5 6.150671123358313e-6]))

calibration(data::Collection{<:UPbSIMSData}, standardages::Collection; baseline::Number=0)

Create a UPbSIMSCalibration object calib given a dataset of standard SIMS analyses data with known ages standardages, to construct a Pb/U RSF vs U/UO2 calibration line.

The resulting calibraiton line is stored in calib.line, while the raw data is stored in calib.data.

Examples

julia> standards = importsimsdata("./examples/data")
10-element Vector{Isoplot.UPbSIMSData{Float64}}:
 Isoplot.UPbSIMSData{Float64}(...)
 Isoplot.UPbSIMSData{Float64}(...)
 ⋮
 Isoplot.UPbSIMSData{Float64}(...)

julia> standardages = fill(1099., length(standards));

julia> calib = calibration(standards, standardages)
Isoplot.UPbSIMSCalibration{Float64}(...)

julia> calib.line
YorkFit{Float64}:
 Least-squares linear fit of the form y = a + bx with
  intercept: 1.25 ± 0.19 (1σ)
  slope    : 0.82 ± 0.066 (1σ)
  MSWD     : 0.2704923359277713

chauvenet(μ::AbstractVector, σ::AbstractVector)

Apply Chauvenet's criterion to a set of data to exclude outliers.

The function calculates the z-scores of the data points, and then calculates the probability p of observing a value as extreme as the z-score under the assumption of normal distribution. It then applies Chauvenet's criterion, marking any data point as an outlier if 2 * N * p < 1.0, where N is the total number of data points.

Examples

dist_ll(dist::Collection, mu::Collection, sigma::Collection, tmin::Number, tmax::Number)
dist_ll(dist::Collection, analyses::Collection{<:Measurement}, tmin::Number, tmax::Number)

Return the log-likelihood of a set of mineral ages with means mu and uncertianty sigma being drawn from a given source (i.e., crystallization / closure) distribution dist, with terms to prevent runaway at low N.

Examples

mu, sigma = collect(100:0.1:101), 0.01*ones(11)
ll = dist_ll(MeltsVolcanicZirconDistribution, mu, sigma, 100, 101)

importsimsdata(dir::String; T=Float64)

Import all .asc files in the directory dir into a vector of UPbSIMSData objects.

Examples

julia> standards = importsimsdata("./examples/data")
10-element Vector{Isoplot.UPbSIMSData{Float64}}:
 Isoplot.UPbSIMSData{Float64}(...)
 Isoplot.UPbSIMSData{Float64}(...)
 ⋮
 Isoplot.UPbSIMSData{Float64}(...)

(a,b) = lsqfit(x::AbstractVector, y::AbstractVector)

Returns the coefficients for a simple linear least-squares regression of the form y = a + bx

Examples

julia> a, b = lsqfit(1:10, 1:10)
2-element Vector{Float64}:
 -1.19542133983862e-15
  1.0

julia> isapprox(a, 0, atol = 1e-12)
true

julia> isapprox(b, 1, atol = 1e-12)
true

metropolis_min!(tmindist::DenseArray, dist::Collection, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
metropolis_min!(tmindist::DenseArray, t0dist::DenseArray, dist::Collection, analyses::Collection{<:UPbAnalysis}; burnin::Integer=0) where {T}

In-place (non-allocating) version of metropolis_min, fills existing array tmindist.

Run a Metropolis sampler to estimate the minimum of a finite-range source distribution dist using samples drawn from that distribution – e.g., estimate zircon eruption ages from a distribution of zircon crystallization ages.

Examples

metropolis_min!(tmindist, 2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)

metropolis_min(nsteps::Integer, dist::Collection, data::Collection{<:Measurement}; burnin::Integer=0, t0prior=Uniform(0,minimum(value.(age68.(analyses)))), lossprior=Uniform(0,100))
metropolis_min(nsteps::Integer, dist::Collection, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
metropolis_min(nsteps::Integer, dist::Collection, analyses::Collection{<:UPbAnalysis; burnin::Integer=0)

Run a Metropolis sampler to estimate the minimum of a finite-range source distribution dist using samples drawn from that distribution – e.g., estimate zircon eruption ages from a distribution of zircon crystallization ages.

Examples

tmindist = metropolis_min(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)

tmindist, t0dist = metropolis_min(2*10^5, HalfNormalDistribution, analyses, burnin=10^5)

metropolis_minmax!(tmindist, tmaxdist, lldist, acceptancedist, nsteps::Integer, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)
metropolis_minmax!(tmindist, tmaxdist, t0dist, lldist, acceptancedist, nsteps::Integer, dist::Collection, analyses::Collection{<:UPbAnalysis}; burnin::Integer=0)

In-place (non-allocating) version of metropolis_minmax, filling existing arrays

Run a Metropolis sampler to estimate the extrema of a finite-range source distribution dist using samples drawn from that distribution – e.g., estimate zircon saturation and eruption ages from a distribution of zircon crystallization ages.

Examples

metropolis_minmax!(tmindist, tmaxdist, lldist, acceptancedist, 2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)

metropolis_minmax(nsteps::Integer, dist::Collection, data::Collection{<:Measurement}; burnin::Integer=0)
metropolis_minmax(nsteps::Integer, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)

Run a Metropolis sampler to estimate the extrema of a finite-range source distribution dist using samples drawn from that distribution – e.g., estimate zircon saturation and eruption ages from a distribution of zircon crystallization ages.

Examples

tmindist, tmaxdist, lldist, acceptancedist = metropolis_minmax(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)

mswd(μ, σ)
mswd(μ ± σ)

Return the Mean Square of Weighted Deviates (AKA the reduced chi-squared statistic) of a dataset with values x and one-sigma uncertainties σ

Examples

julia> x = randn(10)
10-element Vector{Float64}:
 -0.977227094347237
  2.605603343967434
 -0.6869683962845955
 -1.0435377148872693
 -1.0171093080088411
  0.12776158554629713
 -0.7298235147864734
 -0.3164914095249262
 -1.44052961622873
  0.5515207382660242

julia> mswd(x, ones(10))
1.3901517474017941

yorkfit(x::Collection, σx::Collection, y::Collection, σy::Collection, [r])
yorkfit(x::Collection{<:Measurement}, y::Collection{<:Measurement}, [r])
yorkfit(d::Collection{<:AbstractAnalysis})

Uses the York (1968) two-dimensional least-squares fit to calculate a, b, and uncertanties σa, σb for the equation y = a + bx, given x, y, uncertaintes σx, and σy, and optially covarances r.

For further reference, see: York, Derek (1968) "Least squares fitting of a straight line with correlated errors" Earth and Planetary Science Letters 5, 320-324. doi: 10.1016/S0012-821X(68)80059-7

Examples

julia> x = (1:100) .+ randn.();

julia> y = 2*(1:100) .+ randn.();

julia> yorkfit(x, ones(100), y, ones(100))
YorkFit{Float64}:
Least-squares linear fit of the form y = a + bx where
  intercept a : -0.29 ± 0.2 (1σ)
  slope b     : 2.0072 ± 0.0035 (1σ)
  MSWD        : 0.8136665223891004