Export Radiocarbon distributions to Isoplot.jl

Author: brenhinkeller

Project.toml

@@ -20,7 +20,7 @@ StatGeochemBase = "61e559cd-58b4-4257-8528-26bb26ff2b9a"
 [compat]
 DelimitedFiles = "1"
 Distributions = "0.15 - 0.25"
-Isoplot = "0.4"
+Isoplot = "0.5"
 KernelDensity = "0.4 - 0.6"
 LsqFit = "0.7 - 0.15"
 NaNStatistics = "0.6.40"

src/Chron.jl

@@ -39,8 +39,6 @@ module Chron
     include("Utilities.jl")
     # Functions for propagating systematic uncertainties
     include("Systematic.jl")
-    # Intcal2013 calibration curve for radiocarbion
-    include("Intcal.jl")
     # Functions for estimating extrema of a finite-range distribution
     include("DistMetropolis.jl")
     # Functions for stratigraphic modelling
@@ -60,7 +58,7 @@ module Chron
         screen_outliers, BootstrapCrystDistributionKDE
 
     # Custom distributions and related functions
-    export strat_ll, BilinearExponential, Radiocarbon
+    export strat_ll, BilinearExponential
 
     # Dealing with systematic uncertainty
     export add_systematic_uncert_UPb, add_systematic_uncert_ArAr, add_systematic_uncert_UTh, SystematicUncertainty

src/Intcal.jl

@@ -158,110 +158,7 @@
     Base.:+(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A, d.loc + c, d.scl, d.shp, d.skw)
     Base.:*(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A-log(abs(c)), d.loc * c, d.scl*abs(c), d.shp, d.skw)
 
-## --- Radiocarbon distribution type
-
-    struct Radiocarbon{T<:Real} <: ContinuousUnivariateDistribution
-        μ::T
-        σ::T
-        dist::Vector{T}
-        ldist::Vector{T}
-        lcdist::Vector{T}
-        lccdist::Vector{T}
-    end
-
-    function Radiocarbon(μ::T, σ::T, ldist::Vector{T}) where {T<:Real}
-        # Ensure normalization
-        dist = exp.(ldist)
-        normconst = sum(dist) * one(T)
-        dist ./= normconst
-        ldist .-= normconst
-
-        # Cumulative distributions
-        lcdist = nanlogcumsumexp(ldist)
-        lccdist = nanlogcumsumexp(ldist, reverse=true)
 
-        return Radiocarbon{T}(μ, σ, dist, ldist, lcdist, lccdist)
-    end
-
-    function Radiocarbon(Age_14C::Real, Age_14C_sigma::Real, calibration::NamedTuple=intcal13)
-        @assert calibration.Age_Calendar == 1:1:length(calibration.Age_14C)
-        @assert step(calibration.Age_Calendar) == calibration.dt == 1
-
-        ldist = normlogproduct.(Age_14C, Age_14C_sigma, calibration.Age_14C, calibration.Age_sigma)
-
-        # Ensure normalization
-        dist = exp.(ldist)
-        normconst = sum(dist) * calibration.dt
-        dist ./= normconst
-        ldist .-= normconst
-
-        μ = histmean(dist, calibration.Age_Calendar)
-        σ = histstd(dist, calibration.Age_Calendar, corrected=false)
-
-        # Cumulative distributions
-        lcdist = nanlogcumsumexp(ldist)
-        lccdist = nanlogcumsumexp(ldist, reverse=true)
-
-        return Radiocarbon(μ, σ, dist, ldist, lcdist, lccdist)
-    end
-
-    ## Conversions
-    Base.convert(::Type{Radiocarbon{T}}, d::Radiocarbon) where {T<:Real} = Radiocarbon{T}(T(d.μ), T(d.σ), T.(d.ldist))
-    Base.convert(::Type{Radiocarbon{T}}, d::Radiocarbon{T}) where {T<:Real} = d
-
-    ## Parameters
-    Distributions.params(d::Radiocarbon) = (d.dist, d.ldist)
-    @inline Distributions.partype(d::Radiocarbon{T}) where {T<:Real} = T
-
-    Distributions.location(d::Radiocarbon) = d.μ
-    Distributions.scale(d::Radiocarbon) = d.σ
-
-    Base.eltype(::Type{Radiocarbon{T}}) where {T} = T
-
-    ## Evaluation
-    @inline Distributions.pdf(d::Radiocarbon, x::Real) = exp(logpdf(d, x))
-    @inline function Distributions.pdf(d::Radiocarbon{T}, x::Real) where {T}
-        if firstindex(d.ldist) <= x <= lastindex(d.ldist)
-            linterp_at_index(d.dist, x, -maxintfloat(T))
-        else
-            # Treat as a single Normal distrbution if outside of the calibration range
-            pdf(Normal(d.μ, d.σ), x)
-        end
-    end
-    @inline function Distributions.logpdf(d::Radiocarbon{T}, x::Real) where {T}
-        if firstindex(d.ldist) <= x <= lastindex(d.ldist)
-            linterp_at_index(d.ldist, x, -maxintfloat(T))
-        else
-            # Treat as a single Normal distrbution if outside of the calibration range
-            logpdf(Normal(d.μ, d.σ), x)
-        end
-    end
-    @inline Distributions.cdf(d::Radiocarbon, x::Real) = exp(logcdf(d, x))
-    @inline function Distributions.logcdf(d::Radiocarbon{T}, x::Real) where {T}
-        if firstindex(d.lcdist) <= x <= lastindex(d.lcdist)
-            linterp_at_index(d.lcdist, x, -maxintfloat(T))
-        else
-            # Treat as a single Normal distrbution if outside of the calibration range
-            logcdf(Normal(d.μ, d.σ), x)
-        end
-    end
-    @inline Distributions.ccdf(d::Radiocarbon, x::Real) = exp(logccdf(d, x))
-    @inline function Distributions.logccdf(d::Radiocarbon{T}, x::Real) where {T}
-        if firstindex(d.lccdist) <= x <= lastindex(d.lccdist)
-            linterp_at_index(d.lccdist, x, -maxintfloat(T))
-        else
-            # Treat as a single Normal distrbution if outside of the calibration range
-            logccdf(Normal(d.μ, d.σ), x)
-        end
-    end
-
-    ## Statistics
-    Distributions.mean(d::Radiocarbon) = d.μ
-    Distributions.var(d::Radiocarbon) = d.σ * d.σ
-    Distributions.std(d::Radiocarbon) = d.σ
-    Distributions.skewness(d::Radiocarbon) = histskewness(d.dist, eachindex(d.dist), corrected=false)
-    Distributions.kurtosis(d::Radiocarbon) = histkurtosis(d.dist, eachindex(d.dist), corrected=false)
-
 ## --- Extend `normpdf_ll` to deal with Distributions.Normal 
 
     function StatGeochemBase.normpdf_ll(x::AbstractVector, ages::AbstractVector{Normal{T}}) where T

src/Utilities.jl

@@ -26,36 +26,14 @@
     @test skewness(d) ≈ -0.16449904171203025 atol = 1e-6
     @test kurtosis(d) ≈ 0.0865922917156188 atol = 1e-6
 
-## --- Test Radiocarbon distribution
-
-    d = Radiocarbon(1000, 10, intcal13)
-    @test d isa Radiocarbon{Float64}
-    @test pdf(d, 950) ≈ 0.0006472245090905895
-    @test logpdf(d, 950) ≈ -7.342817323513984
-    @test pdf.(d, [900, 950, 1000]) ≈ [6.333126104255167e-8, 0.0006472245090905895, 2.652591331229418e-13]
-    @test cdf.(d, [-1, 900, 950, 1000, 1e5]) ≈ [0.0, 0.000563737957020387, 0.18096988328291203, 0.18312332434946413, 1.0]
-    @test ccdf.(d, [-1, 900, 950, 1000, 1e5]) ≈ [1.0, 0.18255964979458322, 0.0028006656465210025, 7.114355524701459e-11, 0.0]
-    @test logcdf.(d, [-1, 900, 950, 1000, 1e5]) ≈ [-7649.088264850034, -7.480921029645386, -1.7094246522629235, -1.6975954495627632, -0.0] 
-    @test logccdf.(d, [-1, 900, 950, 1000, 1e5]) ≈ [-0.0, -1.700678311173327, -5.8778981591541335, -23.366321375298313, -8.706770436253259e7]
-
-    @test eltype(d) === partype(d) === Float64
-    @test location(d) ≈ 927.2556461305643
-    @test scale(d) ≈ 7.5077650707247
-
-    @test mean(d) ≈ 927.2556461305643
-    @test var(d) ≈ 7.5077650707247^2
-    @test std(d) ≈ 7.5077650707247
-    @test skewness(d) ≈ -4.6327184250788696
-    @test kurtosis(d) ≈ 67.68860469654331
-
 ## --- Other utility functions for log likelihoods
 
     @test strat_ll([66.02,], [BilinearExponential(66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)]) ≈  (-3.251727600957773 -4.265260194845626)
-    @test strat_ll([950,], [Radiocarbon(1000, 10, intcal13)]) ≈ -7.342817323513984
+    @test strat_ll([950,], [Radiocarbon(1000, 10, calibration=intcal13)]) ≈ -7.342817323513984
 
     @test strat_ll([0.0, 0.0], [Normal(0,1), Normal(0,1)]) ≈ 0
     @test strat_ll([0.0, 0.5], [Normal(0,1), Uniform(0,1)]) ≈ 0
     @test strat_ll([0.0, 0.5, 66.02], [Normal(0,1), Uniform(0,1), BilinearExponential(66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)]) ≈ (-3.251727600957773 -4.265260194845626)
-    @test strat_ll([0.0, 0.5, 66.02, 900], [Normal(0,1), Uniform(0,1), BilinearExponential(66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505), Radiocarbon(1000, 10, intcal13)]) ≈ -24.09187457017352
+    @test strat_ll([0.0, 0.5, 66.02, 900], [Normal(0,1), Uniform(0,1), BilinearExponential(66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505), Radiocarbon(1000, 10, calibration=intcal13)]) ≈ -24.09187457017352
 
 ## ---