[Internals] Switch to TaylorDiff.jl (#16)

Author: lrnv

.github/workflows/CI.yml

@@ -19,15 +19,15 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - 'lts'
           - '1'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
       - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

Project.toml

@@ -1,17 +1,17 @@
 name = "WilliamsonTransforms"
 uuid = "48feb556-9bdd-43a2-8e10-96100ec25e22"
+version = "0.2"
 authors = ["Oskar Laverny <[email protected]> and contributors"]
-version = "0.1.7"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
 
 [compat]
 Distributions = "0.25"
 Roots = "1, 2"
-TaylorSeries = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
+TaylorDiff = "0.3"
 julia = "1"
 
 [extras]

docs/Project.toml

@@ -1,3 +1,4 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"

src/WilliamsonTransforms.jl

@@ -7,11 +7,9 @@ module WilliamsonTransforms
 end WilliamsonTransforms
 
 import Distributions
-import TaylorSeries
-import Base.minimum
+import TaylorDiff
 import Roots
 import Base: minimum, maximum
-
 export 𝒲, 𝒲₋₁
 
 """
@@ -36,29 +34,40 @@ References:
 - Williamson, R. E. (1956). Multiply monotone functions and their Laplace transforms. Duke Math. J. 23 189–207. MR0077581
 - McNeil, Alexander J., and Johanna Nešlehová. "Multivariate Archimedean copulas, d-monotone functions and ℓ 1-norm symmetric distributions." (2009): 3059-3097.
 """
-struct 𝒲{TX}
+struct 𝒲{TX, d}
     X::TX
-    d::Int
-    # E::TE
-    function 𝒲(X::TX,d) where TX<:Distributions.UnivariateDistribution
+    function 𝒲(X::TX, ::Val{d}) where {TX<:Distributions.UnivariateDistribution, d}
         @assert minimum(X) ≥ 0 && maximum(X) ≤ Inf 
         @assert d ≥ 2 && isinteger(d) 
-        return new{typeof(X)}(X,d)
+        return new{typeof(X), d}(X)
     end
+    𝒲(X, d::Int) = 𝒲(X, Val(d))
 end
 
-function (ϕ::𝒲)(x)
+function (ϕ::𝒲{TX, d})(x) where {TX,d}
     if x <= 0
         return 1 - Distributions.cdf(ϕ.X,0)
     else
-        return Distributions.expectation(y -> (1 - x/y)^(ϕ.d-1) * (y > x), ϕ.X)
-        # We need to compute the expectation of (1 - x/X)^{d-1}
-        # return ϕ.E(y -> (y > x) * (1 - x/y)^(ϕ.d-1))
+        return Distributions.expectation(y -> (1 - x/y)^(d-1) * (y > x), ϕ.X)
     end
 end
 
-function taylor(f, x, d)
-    return f(x + TaylorSeries.Taylor1([zero(x), one(x)],d)).coeffs
+"""
+    taylor(f::F, x₀, ::Val{d}) where {F,d}
+
+Compute the Taylor series expansion of the function `f` around the point `x₀` up to order `d`, and gives you back all the successive derivatives. 
+
+# Arguments
+- `f`: A function to be expanded.
+- `x₀`: The point around which to expand the Taylor series.
+- `d`: The order up to which the Taylor series is computed.
+
+# Returns
+A tuple with value ``(f(x₀), f'(x₀),...,f^{(d)}(x₀))``.
+"""
+function taylor(f::F, x₀, D::Val{d}) where {F,d} 
+    r = TaylorDiff.derivatives(f, x₀, one(x₀), D)
+    return (r.value, r.partials...)
 end
 
 """
@@ -87,38 +96,36 @@ References:
     - Williamson, R. E. (1956). Multiply monotone functions and their Laplace transforms. Duke Math. J. 23 189–207. MR0077581
     - McNeil, Alexander J., and Johanna Nešlehová. "Multivariate Archimedean copulas, d-monotone functions and ℓ 1-norm symmetric distributions." (2009): 3059-3097.
 """
-struct 𝒲₋₁{Tϕ} <: Distributions.ContinuousUnivariateDistribution
+struct 𝒲₋₁{Tϕ, d} <: Distributions.ContinuousUnivariateDistribution
     ϕ::Tϕ
-    d::Int
-    function 𝒲₋₁(ϕ,d)
+    function 𝒲₋₁(ϕ, ::Val{d}) where d
         @assert ϕ(0.0) == 1.0
         @assert ϕ(float(Inf)) == 0.0
-        # And assertion about d-monotony... how can this be check ? this is hard. 
-        return new{typeof(ϕ)}(ϕ,d)
+        @assert isinteger(d)
+        return new{typeof(ϕ),d}(ϕ)
     end
+    𝒲₋₁(ϕ, d::Int) = 𝒲₋₁(ϕ, Val(d))
 end
-function Distributions.cdf(d::𝒲₋₁, x)
-    rez = zero(x)
-    c_ϕ = taylor(d.ϕ, x, d.d)
-    c_ϕ[end] = max(c_ϕ[end], 0)
-    for k in 0:(d.d-1)
-        if c_ϕ[k+1] != 0 # We need c_ϕ = 0 to dominate x = Inf
-            rez += (-1)^k * x^k * c_ϕ[k+1]
-        end
+function Distributions.cdf(dist::𝒲₋₁{Tϕ, d}, x) where {Tϕ, d}
+    x ≤ 0 && return zero(x)
+    rez, x_pow = zero(x), one(x)
+    c = taylor(dist.ϕ, x, Val(d-1))
+    for k in 1:d
+        rez += iszero(c[k]) ? 0 : x_pow * c[k]
+        x_pow *= -x
     end
-    # simple hack to ensure convergence :
     return isnan(rez) ? one(x) : 1 - rez
-    # return 1-rez
 end
-Distributions.logpdf(d::𝒲₋₁, x::Real) = log(max(0,taylor(x -> Distributions.cdf(d,x), x, 1)[end]))
-_quantile(d::𝒲₋₁, p) = Roots.find_zero(x -> (Distributions.cdf(d, x) - p), (0.0, Inf))
-Distributions.rand(rng::Distributions.AbstractRNG, d::𝒲₋₁) = _quantile(d,rand(rng))
+
+Distributions.logpdf(dist::𝒲₋₁{Tϕ, d}, x) where {Tϕ, d} = log(max(0,taylor(x -> Distributions.cdf(dist,x), x, Val(1))[end]))
+_quantile(dist::𝒲₋₁, p) = Roots.find_zero(x -> (Distributions.cdf(dist, x) - p), (0.0, Inf))
+Distributions.rand(rng::Distributions.AbstractRNG, dist::𝒲₋₁) = _quantile(dist, rand(rng))
 Base.minimum(::𝒲₋₁) = 0.0
 Base.maximum(::𝒲₋₁) = Inf
-function Distributions.quantile(d::𝒲₋₁, p::Real)
-# Validate that p is in the range [0, 1]
+function Distributions.quantile(dist::𝒲₋₁, p::Real)
+    # Validate that p is in the range [0, 1]
     @assert 0 <= p <= 1
-    return _quantile(d,p)
+    return _quantile(dist, p)
 end
 end
 

test/testing_the_paper.jl

@@ -127,5 +127,5 @@ end
 @testitem "testing fractional-dimensional williamson transformation" begin
     using Distributions
     ϕ(x) = exp(-x)
-    @test_throws InexactError 𝒲₋₁(ϕ, 0.7)
+    @test_throws MethodError 𝒲₋₁(ϕ, 0.7)
 end