Fully transition BilinearExpoenential to the Distribution-style system

src/Chron.jl

@@ -55,12 +55,16 @@ module Chron
     export  StratMetropolis, StratMetropolisDist, StratMetropolis14C,
         tMinDistMetropolis, metropolis_min!, metropolis_min,
         metropolis_minmax!, metropolis_minmax,
-        bilinear_exponential, bilinear_exponential_ll,
         screen_outliers, BootstrapCrystDistributionKDE
 
+    # Custom distributions and related functions
+    export strat_ll, BilinearExponential, Radiocarbon
+
+    # Dealing with systematic uncertainty
+    export add_systematic_uncert_UPb, add_systematic_uncert_ArAr, add_systematic_uncert_UTh, SystematicUncertainty
 
-    # Distributions
-    # Mostly already exported by reexporting Isoplot.jl
+    # Relative crystallization/closure distributions
+    # (mostly already exported by reexporting Isoplot.jl)
     export TruncatedNormalDistribution, EllisDistribution, ArClosureDistribution
 
 end # module

src/Fitplot.jl

@@ -137,14 +137,50 @@
 
 ## --- Fit and plot results from stationary distribution of depostion/eruption age distribution model
 
+    # """
+    # ```Julia
+    # bilinear_exponential(x::Number, p::AbstractVector)
+    # ```
+    # Evaluate the value of a "bilinear exponential" function defined by the parameters
+    # `p` at point `x`. This function, which can be used to approximate various asymmetric
+    # probability distributions found in geochronology, is defined by an exponential
+    # function with two log-linear segments joined by an `atan` sigmoid:
+    # ```math
+    # ℯ^{p_1} * ℯ^{v x_s p_4 p_5 - (1-v) x_s p_4/p_5}
+    # ```
+    # where
+    # ```math
+    # v = 1/2 - atan(x_s)/π
+    # ```
+    # is a sigmoid, positive on the left-hand side, and
+    # ```math
+    # x_s = (x - p_2)/p_3
+    # ```
+    # is `x` scaled by mean and standard deviation.
+
+    # The elements of the parameter array `p` may be considered to approximately represent\n
+    #     p[1] # pre-exponential (normaliation constant)
+    #     p[2] # mean (central moment)
+    #     p[3] # standard deviation
+    #     p[4] # sharpness
+    #     p[5] # skew
+
+    # where all parameters `pᵢ` must be nonnegative.
+    # """
+    function bilinear_exponential(x, p)
+        i₀ = firstindex(p)
+        d = BilinearExponential(p[i₀], p[i₀+1], abs(p[i₀+2]), abs(p[i₀+3]), abs(p[i₀+4]))
+        return pdf.(d, x)
+    end
+
     """
     ```julia
     tMinDistMetropolis(smpl::ChronAgeData, nsteps::Int, burnin::Int, dist::Array{Float64})
     ```
     Calculate the minimum limiting (eruption/deposition) age of each sample defined
     in the `smpl` struct, using the `Isoplot.metropolis_min` function, assuming mineral
     ages for each sample are drawn from the source distribution `dist`. Fits a
-    `bilinear_exponential` function to the resulting stationary distribution
+    `BilinearExponential` function to the resulting stationary distribution
     for each sample and stores the results in `smpl.Params` for use by the
     `StratMetropolisDist` function.
 

src/StratMetropolis.jl

@@ -144,7 +144,7 @@
     ```
     Runs the main Chron.jl age-depth model routine for a stratigraphic set of
     samples defined by sample heights and fitted asymmetric age distributions
-    (`bilinear_exponential`) in the `smpl` struct, and an age-depth model
+    (`BilinearExponential`) in the `smpl` struct, and an age-depth model
     configuration defined by the `config` struct.
 
     Optionally, if a `hiatus` struct is provided, the model will additionally
@@ -434,19 +434,19 @@
 ## --- # Internals of the Markov chain
 
     # Use dispatch to let us reduce duplication
-    strat_ll(x, ages::AbstractVector{<:BilinearExponential}) = bilinear_exponential_ll(x, ages)
     strat_ll(x, ages::AbstractVector{<:Radiocarbon}) = interpolate_ll(x, ages)
     strat_ll(x, ages::AbstractVector{<:Normal}) = normpdf_ll(x, ages)
-    function strat_ll(x, ages::AbstractVector)
+    strat_ll(x::Real, age::Distribution) = logpdf(age, x)
+    function strat_ll(x, ages)
         ll = zero(float(eltype(x)))
-        @inbounds for i in eachindex(x,ages)
+        @inbounds for i in eachindex(x, ages)
             ll += logpdf(ages[i], x[i])
         end
         return ll
     end
 
-    adjust!(ages::AbstractVector, chronometer, systematic::Nothing) = ages
-    function adjust!(ages::AbstractVector, chronometer, systematic::SystematicUncertainty) where T
+    adjust!(ages, chronometer, systematic::Nothing) = ages
+    function adjust!(ages, chronometer, systematic::SystematicUncertainty)
         systUPb = randn()*systematic.UPb
         systArAr = randn()*systematic.ArAr
         @assert eachindex(ages)==eachindex(chronometer)

src/Systematic.jl

@@ -15,8 +15,6 @@
         return log.(r .* σtracer .+ 1) ./ λ238_
     end
 
-    export add_systematic_uncert_UPb
-
     function add_systematic_uncert_ArAr(agedistmyr::Vector{<:AbstractFloat}; constants=:Renne)
         if constants==:Min
             # Min et al., 200 compilation
@@ -118,10 +116,6 @@
         return log.(Rdist.*κ_.*λ_./λϵ_ .+ 1)./(λ_.*1000000)
     end
 
-    export add_systematic_uncert_ArAr
-
-## --- End of File
-
     function add_systematic_uncert_UTh(agedistmyr::Vector{<:AbstractFloat})
         # Age = -log(1-Th230ₐ/U238ₐ))/λ230Th
         # Th230ₐ/U238ₐ = 1-exp(-Age*λ230Th)
@@ -160,4 +154,5 @@
         # Maximum possible activity ratio is 1.0, corresponding to secular equilibrium
         return @. -log(1 - min(sTh230_U238_dist*λ230_/λ238_, 1.0)) / λ230_
     end
-    export add_systematic_uncert_UTh
+
+## --- End of File

src/Utilities.jl

@@ -1,5 +1,42 @@
-## --- Fitting non-Gaussian distributions
+## --- Dealing with arbitrary distributions
 
+
+## --- "bilinear exponential" distribution type
+
+    """
+    ```Julia
+    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 
+    of Bayesian eruption age estimation). 
+    
+    This "bilinear exponential" distribution, as the name might suggest, is defined as an 
+    exponential function with two log-linear segments, joined by an arctangent sigmoid:
+    ```math
+    ℯ^{A} * ℯ^{v*xₛ*shp*skw - (1-v)*xₛ*shp/skw}
+    ```
+    where
+    ```math
+    v = 1/2 - atan(xₛ)/π
+    ```
+    is a sigmoid, positive on the left-hand side, and
+    ```math
+    xₛ = (x - μ)/σ
+    ```
+    is `x` scaled by the location parameter `μ` and scale parameter `σ`, 
+    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.
+    """
     struct BilinearExponential{T<:Real} <: ContinuousUnivariateDistribution
         A::T
         μ::T
@@ -31,13 +68,13 @@
     Base.eltype(::Type{BilinearExponential{T}}) where {T} = T
 
     ## Evaluation
-    function Distributions.pdf(d::BilinearExponential, x::Real)
+    @inline function Distributions.pdf(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 exp(d.A + d.shp*d.skw*xs*v - d.shp/d.skw*xs*(1-v))
     end
 
-    function Distributions.logpdf(d::BilinearExponential, x::Real)
+    @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)
@@ -47,102 +84,7 @@
     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)
 
-    """
-    ```Julia
-    bilinear_exponential(x::Number, p::AbstractVector)
-    ```
-    Evaluate the value of a "bilinear exponential" function defined by the parameters
-    `p` at point `x`. This function, which can be used to approximate various asymmetric
-    probability distributions found in geochronology, is defined by an exponential
-    function with two log-linear segments joined by an `atan` sigmoid:
-    ```math
-    ℯ^{p_1} * ℯ^{v x_s p_4 p_5 - (1-v) x_s p_4/p_5}
-    ```
-    where
-    ```math
-    v = 1/2 - atan(x_s)/π
-    ```
-    is a sigmoid, positive on the left-hand side, and
-    ```math
-    x_s = (x - p_2)/p_3
-    ```
-    is `x` scaled by mean and standard deviation.
-
-    The elements of the parameter array `p` may be considered to approximately represent\n
-        p[1] # pre-exponential (normaliation constant)
-        p[2] # mean (central moment)
-        p[3] # standard deviation
-        p[4] # sharpness
-        p[5] # skew
-
-    where all parameters `pᵢ` must be nonnegative.
-    """
-    function bilinear_exponential(x::Number, p::AbstractVector{<:Number})
-        @assert length(p) == 5
-        _, μ, σ, shp, skw = (abs(x) for x in p)
-        xs = (x - μ)/σ # X scaled by mean and variance
-        v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
-        return exp(first(p) + shp*skw*xs*v - shp/skw*xs*(1-v))
-    end
-    function bilinear_exponential(x::AbstractVector, p::AbstractVector{<:Number})
-        @assert length(p) == 5 && firstindex(p) == 1
-        result = Array{float(eltype(x))}(undef,size(x))
-        @inbounds for i ∈ eachindex(x)
-            xs = (x[i] - p[2])/abs(p[3]) # X scaled by mean and variance
-            v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
-            shp = abs(p[4])
-            skw = abs(p[5])
-            result[i] = exp(p[1] + shp*skw*xs*v - shp/skw*xs*(1-v))
-        end
-        return result
-    end
-
-    """
-    ```Julia
-    bilinear_exponential_ll(x, p)
-    ```
-    Return the log likelihood corresponding to a `bilinear_exponential` distribution defined
-    by the parameters `p` evaluate at point `x`.
-
-    If `x` is provided as a number and `p` as a vector, a single evaluation will
-    be returned; if `x` is provided as a vector and `p` as a matrix, the sum
-    over all i for each distribution `p[:,i]` evaluated at `x[i]` will be returned
-    instead.
-
-    See also `bilinear_exponential`
-    """
-    function bilinear_exponential_ll(x::Number, p::AbstractVector{T}) where {T<:Number}
-        @assert length(p) == 5
-        isnan(p[2]) && return zero(T)
-        _, μ, σ, shp, skw = (abs(x) for x in p)
-        xs = (x - μ)/σ # X scaled by mean and variance
-        v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
-        return shp*skw*xs*v - shp/skw*xs*(1-v)
-    end
-    function bilinear_exponential_ll(x::AbstractVector, p::AbstractMatrix{T}) where {T<:Number}
-        ll = zero(T)
-        @inbounds for i ∈ eachindex(x)
-            any(isnan, p[2,i]) && continue
-            xs = (x[i]-p[2,i])/abs(p[3,i]) # X scaled by mean and variance
-            v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
-            shp = abs(p[4,i])
-            skw = abs(p[5,i])
-            ll += shp*skw*xs*v - shp/skw*xs*(1-v)
-        end
-        return lles
-    end
-    function bilinear_exponential_ll(x::AbstractVector, ages::AbstractVector{BilinearExponential{T}}) where T
-        ll = zero(T)
-        @inbounds for i ∈ eachindex(x,ages)
-            pᵢ = ages[i]
-            isnan(pᵢ.μ) && continue
-            xs = (x[i]-pᵢ.μ)/(pᵢ.σ) # X scaled by mean and variance
-            v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
-            shp, skw = pᵢ.shp, pᵢ.skw
-            ll += shp*skw*xs*v - shp/skw*xs*(1-v)
-        end
-        return ll
-    end
+## --- Radiocarbon distribution type
 
     struct Radiocarbon{T}
         μ::T

test/testUtilities.jl

@@ -1,25 +1,25 @@
 ## --- Test bilinear exponential function
 
-    @test bilinear_exponential(66.02, [-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505]) ≈ 0.00010243952455548608
-    @test bilinear_exponential([66.02, 66.08, 66.15], [-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505]) ≈ [0.00010243952455548608, 1.0372066408907184e-13, 1.0569878258022094e-28]
-    @test bilinear_exponential_ll(66.02, [-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505]) ≈ -3.251727600957773
-
-    d = Chron.BilinearExponential([-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505])
-    @test d isa Chron.BilinearExponential{Float64}
+    d = BilinearExponential([-5.9345103368858085, 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 eltype(d) === partype(d) === Float64
     @test params(d) == (-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)
-    @test d + 1 isa Chron.BilinearExponential{Float64}
+    @test d + 1 isa BilinearExponential{Float64}
     @test location(d + 1) == 66.00606672179812 + 1
-    @test d * 2 isa Chron.BilinearExponential{Float64}
+    @test d * 2 isa BilinearExponential{Float64}
     @test location(d * 2) == 66.00606672179812 * 2
     @test scale(d * 2) == 0.17739474265630253 * 2
 
 ## --- Other utility functions for log likelihoods
 
-    @test Chron.strat_ll([0.0, 0.0], [Normal(0,1), Normal(0,1)]) ≈ 0
-    @test Chron.strat_ll([0.0, 0.5], [Normal(0,1), Uniform(0,1)]) ≈ -0.9189385332046728
+    @test strat_ll(66.02, BilinearExponential(-5.9345103368858085, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)) ≈  (-3.251727600957773 -5.9345103368858085)
+
+    @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)
 
 ## ---