[Internals] Back to TaylorSeries.jl (#313)

Author: lrnv

Project.toml

@@ -24,7 +24,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
+TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -59,9 +59,9 @@ StableRNGs = "1"
 Statistics = "1"
 StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1.3"
+TaylorSeries = "0.20"
 Test = "1"
 TestItemRunner = "1"
-WilliamsonTransforms = "0.2"
 julia = "1"
 
 [extras]

README.md

@@ -11,7 +11,6 @@
      <a href="https://github.com/lrnv/Copulas.jl/actions/workflows/CI.yml?query=branch%3Amain"><img src="https://github.com/lrnv/Copulas.jl/actions/workflows/CI.yml/badge.svg?branch=main" alt="Build Status" /></a>
      <a href="https://codecov.io/gh/lrnv/Copulas.jl"><img src="https://codecov.io/gh/lrnv/Copulas.jl/branch/main/graph/badge.svg"/></a>
      <a href="https://github.com/JuliaTesting/Aqua.jl"><img src="https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg" alt="Aqua QA" /></a>
-    <!-- <a href="https://benchmark.tansongchen.com/TaylorDiff.jl"><img src="https://img.shields.io/buildkite/2c801728055463e7c8baeeb3cc187b964587235a49b3ed39ab/main.svg?label=benchmark" alt="Benchmark Status" /></a> -->
 <br />
     <a href="https://opensource.org/licenses/MIT"><img src="https://img.shields.io/badge/license-MIT-blue.svg" alt="License: MIT" /></a>
     <a href="https://github.com/SciML/ColPrac"><img src="https://img.shields.io/badge/contributor's%20guide-ColPrac-blueviolet" alt="ColPrac: Contributor's Guide on Collaborative Practices for Community Packages" /></a>

docs/src/bestiary/archimedean.md

@@ -96,7 +96,7 @@ In this package, we implemented it through the [`WilliamsonGenerator`](@ref) cla
 
 `WilliamsonGenerator(X::UnivariateRandomVariable, d)`.
 
-This function computes the Williamson d-transform of the provided random variable $X$ using the [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package. See [williamson1955multiply, mcneil2009](@cite) for the literature. 
+This function computes the Williamson d-transform of the provided random variable $X$. See [williamson1955multiply, mcneil2009](@cite) for the literature. 
 
 !!! info "`max_monotony` of Williamson generators"
     The $d$-transform of a positive random variable is $d$-monotone but not $k$-monotone for any $k > d$. Its max monotony is therefore $d$. This has a few implications, one of the biggest is that the $d$-variate Archimedean copula that corresponds has no density.
@@ -115,15 +115,15 @@ The Williamson d-transform is a bijective transformation[^1] from the set of pos
 
     This bijection is to be taken carefuly: the bijection is between random variables *with unit scales* and generators *with common value at 1*, sicne on both rescaling does not change the underlying copula. 
 
-This transformation is implemented through one method in the Generator interface that is worth talking a bit about : `williamson_dist(G::Generator, d)`. This function computes the inverse Williamson d-transform of the d-monotone archimedean generator ϕ, still using the [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package. See [williamson1955multiply, mcneil2009](@cite).
+This transformation is implemented through one method in the Generator interface that is worth talking a bit about : `williamson_dist(G::Generator, d)`. This function computes the inverse Williamson d-transform of the d-monotone archimedean generator ϕ. See [williamson1955multiply, mcneil2009](@cite).
 
 To put it in a nutshell, for ``\phi`` a ``d``-monotone archimedean generator, the inverse Williamson-d-transform of ``\\phi`` is the cumulative distribution function ``F`` of a non-negative random variable ``R``, defined by : 
 
 ```math
 F(x) = 𝒲_{d}^{-1}(\phi)(x) = 1 - \frac{(-x)^{d-1} \phi_+^{(d-1)}(x)}{k!} - \sum_{k=0}^{d-2} \frac{(-x)^k \phi^{(k)}(x)}{k!}
 ```
 
-The [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package implements this transformation (and its inverse, the Williamson d-transfrom) in all generality. It returns this cumulative distribution function in the form of the corresponding random variable `<:Distributions.ContinuousUnivariateDistribution` from `Distributions.jl`. You may then compute : 
+It returns this cumulative distribution function in the form of the corresponding random variable `<:Distributions.ContinuousUnivariateDistribution` from `Distributions.jl`. You may then compute : 
 * The cdf via `Distributions.cdf`
 * The pdf via `Distributions.pdf` and the logpdf via `Distributions.logpdf`
 * Samples from the distribution via `rand(X,n)`.
@@ -184,14 +184,14 @@ This is why `williamson_dist(G::Generator,d)` is such an important function in t
 
 ```@example
 using Copulas: williamson_dist, FrankGenerator
-williamson_dist(FrankGenerator(7), Val{3}())
+williamson_dist(FrankGenerator(7), 3)
 ```
 
 For the Frank Copula, as for many classic copulas, the distribution used is known. We pull some of them from `Distributions.jl` but implement a few more, as this Logarithmic one. Another useful example are negatively-dependent Clayton copulas: 
 
 ```@example
 using Copulas: williamson_dist, ClaytonGenerator
-williamson_dist(ClaytonGenerator(-0.2), Val{3}())
+williamson_dist(ClaytonGenerator(-0.2), 3)
 ```
 
 for which the corresponding distribution is known but has no particular name, thus we implemented it under the `ClaytonWilliamsonDistribution` name.
@@ -209,7 +209,7 @@ for which the corresponding distribution is known but has no particular name, th
     We use this fraily approach for several generators, since sometimes it is faster, including e.g. the Clayton one with positive dependence:
     ```@example
     using Copulas: williamson_dist, ClaytonGenerator
-    williamson_dist(ClaytonGenerator(10), Val{3}())
+    williamson_dist(ClaytonGenerator(10), 3)
     ```
 
 

docs/src/bestiary/empirical.md

@@ -216,7 +216,7 @@ Usage:
 
 - Build from data `u::d×n` (raw or pseudos): `Ĝ = EmpiricalGenerator(u; pseudo_values=true)`
 - Use directly in an Archimedean copula: `Ĉ = ArchimedeanCopula(d, Ĝ)`
-- Access the fitted radial law: `R̂ = williamson_dist(Ĝ, Val{d}())`
+- Access the fitted radial law: `R̂ = williamson_dist(Ĝ, d)`
 
 ```@docs; canonical=false
 EmpiricalGenerator

docs/src/examples/archimedean_radial_estimation.md

@@ -205,7 +205,7 @@ R = Dirac(1.0)
 d, n = 3, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -216,7 +216,7 @@ R = DiscreteNonParametric([1.0, 4.0, 8.0], fill(1/3, 3))
 d, n = 2, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -227,7 +227,7 @@ R = DiscreteNonParametric([1.0, 4.0, 8.0], fill(1/3, 3))
 d, n = 3, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -238,7 +238,7 @@ R = LogNormal(1, 3)
 d, n = 10, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R, logged=true)
 ```
 
@@ -249,7 +249,7 @@ R = Pareto(1.0, 1/2)
 d, n = 10, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R, logged=true)
 ```
 

src/ArchimaxCopula.jl

@@ -201,7 +201,7 @@ function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
 
     s    = S * A0
     φp   = ϕ⁽¹⁾(C.gen, s)            # < 0
-    φpp  = ϕ⁽ᵏ⁾(C.gen, Val(2), s)            # > 0
+    φpp  = ϕ⁽ᵏ⁾(C.gen, 2, s)            # > 0
 
     base = su*sv + (φp/φpp)*suv
     base > 0 || return T(-Inf)

src/ArchimedeanCopula.jl

@@ -65,12 +65,12 @@ function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG}
     if !all(0 .< u .< 1)
         return eltype(u)(-Inf)
     end
-    return log(max(ϕ⁽ᵏ⁾(C.G, Val{d}(), sum(ϕ⁻¹.(C.G, u))) * prod(ϕ⁻¹⁽¹⁾.(C.G, u)), 0))
+    return log(max(ϕ⁽ᵏ⁾(C.G, d, sum(ϕ⁻¹.(C.G, u))) * prod(ϕ⁻¹⁽¹⁾.(C.G, u)), 0))
 end
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimedeanCopula{d, TG}, x::AbstractVector{T}) where {T<:Real, d, TG}
     # By default, we use the Williamson sampling.
     Random.randexp!(rng,x)
-    r = rand(rng, williamson_dist(C.G, Val{d}()))
+    r = rand(rng, williamson_dist(C.G, d))
     sx = sum(x)
     for i in 1:length(C)
         x[i] = ϕ(C.G,r * x[i]/sx)
@@ -135,7 +135,7 @@ function rosenblatt(C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real}) where
                 if iszero(rⱼ)
                      U[j,i] = zero(rⱼ)
                 else
-                    A, B = ϕ⁽ᵏ⁾(C.G, Val(j - 1), rⱼ), ϕ⁽ᵏ⁾(C.G, Val(j - 1), rⱼ₋₁)
+                    A, B = ϕ⁽ᵏ⁾(C.G, j - 1, rⱼ), ϕ⁽ᵏ⁾(C.G, j - 1, rⱼ₋₁)
                     U[j,i] = A / B
                 end
             end
@@ -155,8 +155,8 @@ function inverse_rosenblatt(C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real
             elseif !isfinite(Cᵢⱼ)
                 U[j,i] = zero(Cᵢⱼ)
             else
-                Dᵢⱼ = ϕ⁽ᵏ⁾(C.G, Val{j - 1}(), Cᵢⱼ) * u[j,i]
-                R = ϕ⁽ᵏ⁾⁻¹(C.G, Val{j - 1}(), Dᵢⱼ; start_at=Cᵢⱼ)
+                Dᵢⱼ = ϕ⁽ᵏ⁾(C.G, j - 1, Cᵢⱼ) * u[j,i]
+                R = ϕ⁽ᵏ⁾⁻¹(C.G, j - 1, Dᵢⱼ; start_at=Cᵢⱼ)
                 U[j, i] = ϕ(C.G, R - Cᵢⱼ)
                 Cᵢⱼ = R
             end
@@ -172,10 +172,10 @@ function DistortionFromCop(C::ArchimedeanCopula, js::NTuple{p,Int}, uⱼₛ::NTu
     @inbounds for u in uⱼₛ
         sJ += ϕ⁻¹(C.G, u)
     end
-    return ArchimedeanDistortion(C.G, p, float(sJ), float(T(ϕ⁽ᵏ⁾(C.G, Val{p}(), sJ))))
+    return ArchimedeanDistortion(C.G, p, float(sJ), float(T(ϕ⁽ᵏ⁾(C.G, p, sJ))))
 end
 function ConditionalCopula(C::ArchimedeanCopula{D}, ::NTuple{p,Int}, uⱼₛ::NTuple{p,Float64}) where {D, p}
-    return ArchimedeanCopula(D - p, TiltedGenerator(C.G, Val{p}(), sum(ϕ⁻¹.(C.G, uⱼₛ))))
+    return ArchimedeanCopula(D - p, TiltedGenerator(C.G, p, sum(ϕ⁻¹.(C.G, uⱼₛ))))
 end
 SubsetCopula(C::ArchimedeanCopula{d,TG}, dims::NTuple{p, Int}) where {d,TG,p} = ArchimedeanCopula(length(dims), C.G)
 

src/Copulas.jl

@@ -12,7 +12,6 @@ module Copulas
     import ForwardDiff
     import HCubature
     import MvNormalCDF
-    import WilliamsonTransforms
     import Combinatorics
     import LogExpFunctions
     import QuadGK
@@ -22,6 +21,7 @@ module Copulas
     import LambertW
     import Optim
     import Printf
+    import TaylorSeries
 
     # Main code
     include("utils.jl")
@@ -82,6 +82,7 @@ module Copulas
     include("EllipticalCopulas/TCopula.jl")
 
     # Archimedean copulas
+    include("WilliamsonTransforms.jl")
     include("Generator.jl")
     include("ArchimedeanCopula.jl")