Add skewness and kurtosis for Radiocarbon, use effective kurtosis

Author: brenhinkeller

Project.toml

@@ -23,7 +23,7 @@ Distributions = "0.15 - 0.25"
 Isoplot = "0.3"
 KernelDensity = "0.4 - 0.6"
 LsqFit = "0.7 - 0.13"
-NaNStatistics = "0.6.39"
+NaNStatistics = "0.6.40"
 Plots = "1"
 ProgressMeter = "1"
 QuadGK = "2"

src/Utilities.jl

@@ -141,7 +141,7 @@
     Distributions.var(d::BilinearExponential; mean=mean(d)) = first(quadgk(x->(x-mean)^2*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))
     Distributions.std(d::BilinearExponential; mean=mean(d)) = sqrt(var(d; mean))
     Distributions.skewness(d::BilinearExponential; mean=mean(d), std=std(d; mean)) = first(quadgk(x->(x-mean)^3*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))/std^3
-    Distributions.kurtosis(d::BilinearExponential; mean=mean(d), std=std(d; mean)) = first(quadgk(x->(x-mean)^4*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))/std^4
+    Distributions.kurtosis(d::BilinearExponential; mean=mean(d), std=std(d; mean)) = first(quadgk(x->(x-mean)^4*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))/std^4 - 3
 
     ## Affine transformations
     Base.:+(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A, d.μ + c, d.σ, d.shp, d.skw)
@@ -152,42 +152,46 @@
     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
-        normconst = sum(exp.(ldist)) * one(T)
+        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}(μ, σ, ldist, lcdist, lccdist)
+        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)
-        dist = exp.(ldist)
 
-        # Normalize
+        # 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(μ, σ, ldist, lcdist, lccdist)
+        return Radiocarbon(μ, σ, dist, ldist, lcdist, lccdist)
     end
 
     ## Conversions
@@ -205,6 +209,14 @@
 
     ## 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))
@@ -236,6 +248,8 @@
     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)
 
 ## --- 
 

test/testUtilities.jl

@@ -24,7 +24,7 @@
     @test var(d) ≈ 0.02080666596420308^2 atol = 1e-6
     @test std(d) ≈ 0.02080666596420308 atol = 1e-6
     @test skewness(d) ≈ -0.16449904171203025 atol = 1e-6
-    @test kurtosis(d) ≈ 3.0865922917156188 atol = 1e-6
+    @test kurtosis(d) ≈ 0.0865922917156188 atol = 1e-6
 
 ## --- Test Radiocarbon distribution
 
@@ -45,6 +45,8 @@
     @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