Isoplot

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

μ = [r²⁰⁷Pb²³⁵U, r²⁰⁶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])

Apply Chauvenet's criterion to a set of data to identify 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.

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)

(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

wμ, wσ, mswd = wmean(μ, σ; corrected=true, chauvenet=false)
wμ ± wσ, mswd = wmean(μ ± σ; corrected=true, chauvenet=false)

The weighted mean, with or without the "geochronologist's MSWD correction" to uncertainty. You may specify your means and standard deviations either as separate vectors μ and σ, or as a single vector x of Measurements equivalent to x = μ .± σ

In all cases, σ is assumed to reported as actual sigma (i.e., 1-sigma).

If corrected=true, the resulting uncertainty of the weighted mean is expanded by a factor of sqrt(mswd) to attempt to account for dispersion dispersion when the MSWD is greater than 1

If chauvenet=true, outliers will be removed before the computation of the weighted mean using Chauvenet's criterion.

Examples

julia> x = randn(10)
10-element Vector{Float64}:
  0.4612989881720301
 -0.7255529837975242
 -0.18473979056481055
 -0.4176427262202118
 -0.21975911391551833
 -1.6250003193791873
 -1.6185557291787287
  0.25315988825847513
 -0.4979804844182867
  1.3565281078086726

julia> y = ones(10);

julia> wmean(x, y)
(-0.321824416323509, 0.31622776601683794, 0.8192171477885678)

julia> wmean(x .± y)
(-0.32 ± 0.32, 0.8192171477885678)

julia> wmean(x .± y./10)
(-0.322 ± 0.032, 81.9217147788568)

julia> wmean(x .± y./10, corrected=true)
(-0.32 ± 0.29, 81.9217147788568)

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