Add cdf,ccdf and other functions for BilinearExponential, using…

… QuadGK.jl

Project.toml

@@ -1,7 +1,7 @@
 name = "Chron"
 uuid = "68885b1f-77b5-52a7-b2e7-6a8014c56b98"
 authors = ["C. Brenhin Keller <cbkeller@dartmouth.edu>"]
-version = "0.5.4"
+version = "0.5.5"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -12,6 +12,7 @@ LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 NaNStatistics = "b946abbf-3ea7-4610-9019-9858bfdeaf2d"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatGeochemBase = "61e559cd-58b4-4257-8528-26bb26ff2b9a"
@@ -25,6 +26,7 @@ LsqFit = "0.7 - 0.13"
 NaNStatistics = "0.6.38"
 Plots = "1"
 ProgressMeter = "1"
+QuadGK = "2"
 Reexport = "0.2, 1.0"
 SpecialFunctions = "0.10, 1, 2"
 StatGeochemBase = "0.5.9, 0.6"

src/Chron.jl

@@ -32,6 +32,7 @@ module Chron
     using KernelDensity: kde
     using DelimitedFiles
     using Distributions
+    using QuadGK
 
     include("Utilities.jl")
     # Functions for propagating systematic uncertainties

src/Fitplot.jl

@@ -272,7 +272,10 @@
 
                     # Fit nonlinear model
                     fobj = curve_fit(bilinear_exponential,bincenters,N,p)
-                    fobj.param[2:end] .= abs.(fobj.param[2:end]) # Ensure positive
+                    fobj.param[3:end] .= abs.(fobj.param[3:end]) # Ensure positive
+                    μ, σ, shp = fobj.param[2:4]
+                    area, e = quadgk(x->bilinear_exponential(x,fobj.param), μ-200σ/shp, μ+200σ/shp)
+                    fobj.param[1] -= log(area) # Normalize
                     smpl.Params[:,i] = fobj.param
                 end
 
@@ -332,7 +335,9 @@
                 # Fit nonlinear model
                 x = μ .+ (-10σ:σ/10:10σ)
                 fobj = curve_fit(bilinear_exponential, x, normpdf(μ, σ, x), p)
-                fobj.param[2:end] .= abs.(fobj.param[2:end]) # Ensure positive
+                area, e = quadgk(x->bilinear_exponential(x,fobj.param), μ-200σ, μ+200σ)
+                fobj.param[1] -= log(area) # Normalize
+                fobj.param[3:end] .= abs.(fobj.param[3:end]) # Ensure positive
                 fobj.param[5] = 1 # Must be symmetrical
                 smpl.Params[:,i] = fobj.param
             end

src/Utilities.jl

@@ -5,15 +5,15 @@
 
     """
     ```Julia
+    BilinearExponential(μ, σ, shp, skw)
+    BilinearExponential(p::AbstractVector)
     struct BilinearExponential{T<:Real} <: ContinuousUnivariateDistribution
         A::T
         μ::T
         σ::T
         shp::T
         skw::T
     end
-    BilinearExponential(p::AbstractVector)
-    BilinearExponential(p::NTuple{5,T})
     ```
     A five-parameter pseudo-distribution which can be used to approximate various 
     asymmetric probability distributions found in geochronology (including as a result 
@@ -36,22 +36,33 @@
     In addition to the scale parameter `σ`, the additional shape parameters `shp` and `skw` 
     (which control the sharpness and skew of the resulting distribution, respectively), are 
     all three required to be nonnegative.
+
+    If only four parameters `(μ, σ, shp, skw)` are specified, the normalization constant `A` 
+    will be calculated such that the resulting distribution is normalized.
     """
     struct BilinearExponential{T<:Real} <: ContinuousUnivariateDistribution
         A::T
         μ::T
         σ::T
         shp::T
         skw::T
+        function BilinearExponential(A::T, μ::T, σ::T, shp::T, skw::T) where {T}
+            new{T}(A, μ, abs(σ), abs(shp), abs(skw))
+        end
+        function BilinearExponential{T}(A, μ, σ, shp, skw) where {T<:Real}
+            new{T}(A, μ, abs(σ), abs(shp), abs(skw))
+        end
     end
     function BilinearExponential(p::AbstractVector{T}) where {T}
-        @assert length(p) == 5
-        @assert all(x->!(x<0), Iterators.drop(p,1))
-        BilinearExponential{T}(p...)
+        i₀ = firstindex(p)
+        @assert eachindex(p) == i₀:i₀+4
+        BilinearExponential{T}(p[i₀], p[i₀+1], p[i₀+2], p[i₀+3], p[i₀+4])
     end
-    function BilinearExponential(p::NTuple{5,T}) where {T}
-        @assert all(x->!(x<0), Iterators.drop(p,1))
-        BilinearExponential{T}(p...)
+    function BilinearExponential(μ::T, σ::T, shp::T, skw::T) where {T<:Real}
+        # Calculate A such that the resulting distribution is properly normalized
+        d₀ = BilinearExponential{T}(zero(T), μ, σ, shp, skw)
+        area, e = quadgk(x->pdf(d₀, x), μ-200abs(σ/shp), μ+200abs(σ/shp))
+        BilinearExponential{float(T)}(-log(area), μ, σ, shp, skw)
     end
 
     ## Conversions
@@ -73,16 +84,22 @@
         v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
         return exp(d.A + d.shp*d.skw*xs*v - d.shp/d.skw*xs*(1-v))
     end
-
+    Distributions.cdf(d::BilinearExponential, x::Real) = first(quadgk(x->pdf(d,x), d.μ-200*d.σ/d.shp, x))
+    Distributions.ccdf(d::BilinearExponential, x::Real) = first(quadgk(x->pdf(d,x), x, d.μ+200*d.σ/d.shp))
     @inline function Distributions.logpdf(d::BilinearExponential, x::Real)
         xs = (x - d.μ)/d.σ # X scaled by mean and variance
         v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
         return d.A + d.shp*d.skw*xs*v - d.shp/d.skw*xs*(1-v)
     end
 
+    ## Statistics
+    Distributions.mean(d::BilinearExponential) = first(quadgk(x->x*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))
+    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))
+
     ## Affine transformations
     Base.:+(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A, d.μ + c, d.σ, d.shp, d.skw)
-    Base.:*(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A, d.μ * c, d.σ * abs(c), d.shp, d.skw)
+    Base.:*(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A-log(abs(c)), d.μ * c, d.σ*abs(c), d.shp, d.skw)
 
 ## --- Radiocarbon distribution type
 
@@ -130,13 +147,13 @@
     @inline function Distributions.pdf(d::Radiocarbon{T}, x::Real) where {T}
         return linterp_at_index(d.dist, x, zero(T))
     end
-
     @inline function Distributions.logpdf(d::Radiocarbon{T}, x::Real) where {T}
         return linterp_at_index(d.ldist, x, -maxintfloat(T))
     end
 
     ## Statistics
     Distributions.mean(d::Radiocarbon) = d.μ
+    Distributions.var(d::Radiocarbon) = d.σ * d.σ
     Distributions.std(d::Radiocarbon) = d.σ
 
 ## --- 

test/testUtilities.jl

@@ -1,19 +1,25 @@
 ## --- Test bilinear exponential function
 
-    d = BilinearExponential([-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505])
+    d = BilinearExponential([-4.265260194845627, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505])
     @test d isa BilinearExponential{Float64}
-    @test pdf(d, 66.02) ≈ 0.00010243952455548608
-    @test logpdf(d, 66.02) ≈ (-3.251727600957773 -5.9345103368858085)
-    @test pdf.(d, [66.02, 66.08, 66.15]) ≈ [0.00010243952455548608, 1.0372066408907184e-13, 1.0569878258022094e-28]
+    @test pdf(d, 66.02) ≈ 0.0005437680417927236
+    @test logpdf(d, 66.02) ≈ (-3.251727600957773 -4.265260194845627)
+    @test pdf.(d, [66.02, 66.08, 66.15]) ≈ [0.0005437680417927236, 5.505685686251538e-13, 5.610687893459488e-28]
+    @test cdf.(d, [66.02, 66.08, 66.15]) ≈ [0.9999979646869678, 0.9999999999999903, 1.0000000000000446]
+    @test ccdf.(d, [66.02, 66.08, 66.15]) ≈ [2.035313044102181e-6, 1.2768666553428429e-15, 1.0225965286520566e-30]
 
     @test eltype(d) === partype(d) === Float64
-    @test params(d) == (-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)
+    @test params(d) == (-4.265260194845627, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)
     @test d + 1 isa BilinearExponential{Float64}
     @test location(d + 1) == 66.00606672179812 + 1
     @test d * 2 isa BilinearExponential{Float64}
     @test location(d * 2) == 66.00606672179812 * 2
     @test scale(d * 2) == 0.17739474265630253 * 2
 
+    d = BilinearExponential(66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)
+    @test mean(d) ≈ 65.93219453443876 atol = 1e-6
+    @test var(d) ≈ 0.02080666596420308^2 atol = 1e-6
+    @test std(d) ≈ 0.02080666596420308 atol = 1e-6
 
 ## --- Test Radiocarbon distribution
 
@@ -28,6 +34,7 @@
     @test scale(d) ≈ 7.412902965559571
 
     @test mean(d) ≈ 927.2875322569857
+    @test var(d) ≈ 7.412902965559571^2
     @test std(d) ≈ 7.412902965559571
 
 ## --- Other utility functions for log likelihoods