Improve implementation of Radiocarbon, add efficient logcdf, logccdf, etc.

Author: brenhinkeller

Project.toml

@@ -23,11 +23,11 @@ Distributions = "0.15 - 0.25"
 Isoplot = "0.3"
 KernelDensity = "0.4 - 0.6"
 LsqFit = "0.7 - 0.13"
-NaNStatistics = "0.6.38"
+NaNStatistics = "0.6.39"
 Plots = "1"
 ProgressMeter = "1"
 QuadGK = "2"
 Reexport = "0.2, 1.0"
 SpecialFunctions = "0.10, 1, 2"
-StatGeochemBase = "0.5.9, 0.6"
+StatGeochemBase = "0.6.2"
 julia = "1"

src/Utilities.jl

@@ -152,28 +152,42 @@
     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}
-        dist = exp.(ldist)
-        dist ./= sum(dist) * 1 # Normalize
-        return Radiocarbon{T}(μ, σ, dist, ldist)
+        # Ensure normalization
+        normconst = sum(exp.(ldist)) * one(T)
+        ldist .-= normconst
+
+        # Cumulative distributions
+        lcdist = nanlogcumsumexp(ldist)
+        lccdist = nanlogcumsumexp(ldist, reverse=true)
+    
+        return Radiocarbon{T}(μ, σ, 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 = normproduct_ll.(Age_14C, Age_14C_sigma, calibration.Age_14C, calibration.Age_sigma)
-
+        ldist = normlogproduct.(Age_14C, Age_14C_sigma, calibration.Age_14C, calibration.Age_sigma)
         dist = exp.(ldist)
-        dist ./= sum(dist) * calibration.dt # Normalize
+
+        # Normalize
+        normconst = sum(dist) * calibration.dt
+        dist ./= normconst
+        ldist .-= normconst
         μ = histmean(dist, calibration.Age_Calendar)
         σ = histstd(dist, calibration.Age_Calendar, corrected=false)
 
-        return Radiocarbon(μ, σ, dist, ldist)
+        # Cumulative distributions
+        lcdist = nanlogcumsumexp(ldist)
+        lccdist = nanlogcumsumexp(ldist, reverse=true)
+
+        return Radiocarbon(μ, σ, ldist, lcdist, lccdist)
     end
 
     ## Conversions
@@ -190,11 +204,32 @@
     Base.eltype(::Type{Radiocarbon{T}}) where {T} = T
 
     ## Evaluation
-    @inline function Distributions.pdf(d::Radiocarbon{T}, x::Real) where {T}
-        return linterp_at_index(d.dist, x, zero(T))
-    end
+    @inline Distributions.pdf(d::Radiocarbon, x::Real) = exp(logpdf(d, x))
     @inline function Distributions.logpdf(d::Radiocarbon{T}, x::Real) where {T}
-        return linterp_at_index(d.ldist, x, -maxintfloat(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

test/testRadiocarbon.jl

@@ -37,7 +37,7 @@ for i = 1:nSamples
     smpl.Age_sigma[i] = std(samples)
 
     # Populate smpl.Params with log likelihood for each sample
-    smpl.Params[:,i] = normproduct_ll.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], calibration.Age_14C, calibration.Age_sigma)
+    smpl.Params[:,i] = normlogproduct.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], calibration.Age_14C, calibration.Age_sigma)
 end
 
 # Run stratigraphic model
@@ -111,7 +111,7 @@ for i = 1:nSamples
     smpl.Age_sigma[i] = std(samples)
 
     # Populate smpl.Params with log likelihood for each sample
-    smpl.Params[:,i] = normproduct_ll.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], calibration.Age_14C, calibration.Age_sigma)
+    smpl.Params[:,i] = normlogproduct.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], calibration.Age_14C, calibration.Age_sigma)
 end
 
 # Run stratigraphic model

test/testUtilities.jl

@@ -30,26 +30,30 @@
 
     d = Radiocarbon(1000, 10, intcal13)
     @test d isa Radiocarbon{Float64}
-    @test pdf(d, 950) ≈ 0.0036210273644940918
-    @test logpdf(d, 950) ≈ -3.4446721311475414 
-    @test pdf.(d, [900, 950, 1000]) ≈ [3.372067127114498e-7, 0.0036210273644940918, 1.5582187464154278e-12]
+    @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.2875322569857
-    @test scale(d) ≈ 7.412902965559571
+    @test location(d) ≈ 927.2556461305643
+    @test scale(d) ≈ 7.5077650707247
 
-    @test mean(d) ≈ 927.2875322569857
-    @test var(d) ≈ 7.412902965559571^2
-    @test std(d) ≈ 7.412902965559571
+    @test mean(d) ≈ 927.2556461305643
+    @test var(d) ≈ 7.5077650707247^2
+    @test std(d) ≈ 7.5077650707247
 
 ## --- Other utility functions for log likelihoods
 
     @test strat_ll(66.02, BilinearExponential(-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)) ≈  (-3.251727600957773 -5.9345103368858085)
-    @test strat_ll(950, Radiocarbon(1000, 10, intcal13)) ≈ -3.4446721311475414 
+    @test strat_ll(950, Radiocarbon(1000, 10, 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.9189385332046728
     @test strat_ll([0.0, 0.5, 66.02], [Normal(0,1), Uniform(0,1), BilinearExponential(-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)]) ≈ (-0.9189385332046728 -3.251727600957773 -5.9345103368858085)
-    @test strat_ll([0.0, 0.5, 66.02, 900], [Normal(0,1), Uniform(0,1), BilinearExponential(-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505), Radiocarbon(1000, 10, intcal13)]) ≈ -22.831420814939655
+    @test strat_ll([0.0, 0.5, 66.02, 900], [Normal(0,1), Uniform(0,1), BilinearExponential(-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505), Radiocarbon(1000, 10, intcal13)]) ≈ -26.680063245418374
 
 ## ---