Rewrite -- better version

Project.toml

@@ -5,19 +5,16 @@ version = "0.0.1"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-Expectations = "2fe49d83-0758-5602-8f54-1f90ad0d522b"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
 julia = "1"
 
 [extras]
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "TestItemRunner"]
+test = ["Test", "TestItemRunner", "SpecialFunctions"]

src/WilliamsonTransform.jl

@@ -1,13 +1,9 @@
 module WilliamsonTransform
 
 import Distributions
-import Expectations
 import TaylorSeries
 import Base.minimum
-import SpecialFunctions
 import Roots
-import StatsBase
-
 
 export 𝒲, 𝒲₋₁
 
@@ -28,25 +24,24 @@ For a univariate non-negative random variable ``X`` for distribution function ``
 𝒲_{X,d}(x) = \\int_{x}^{\\infty} \\left(1 - \\frac{x}{t}\\right)^{d-1} dF(t) = \\mathbb E\\left( (1 - \\frac{x}{X})^{d-1}_+\\right) \\mathbb 1_{x > 0} + \\left(1 - F(0)\\right)\\mathbb 1_{x <0}
 ```
 """
-struct 𝒲{TX,TE}
+struct 𝒲{TX}
     X::TX
     d::Int64
-    E::TE
+    # E::TE
     function 𝒲(X::TX,d) where TX<:Distributions.UnivariateDistribution
-        # S = support(X)
-        # @assert S.lb ≥ 0 && S.ub ≤ Inf # check that X is indeed non-negative. 
-        # @assert d ≥ 2 && isinteger(d) # check that d is an integer greater than 2.
-        E = Expectations.expectation(X) 
-        return new{typeof(X),typeof(E)}(X,d,E)
+        @assert minimum(X) ≥ 0 && maximum(X) ≤ Inf 
+        @assert d ≥ 2 && isinteger(d) 
+        return new{typeof(X)}(X,d)
     end
 end
 
 function (ϕ::𝒲)(x)
     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 ϕ.E(y -> (y > x) * (1 - x/y)^(ϕ.d-1))
     end
 end
 
@@ -69,28 +64,32 @@ The cumulative distribution function of this random variable is given by:
 𝒲₋₁(X,d)(x) = 1 - \\frac{(-x)^{d-1} \\phi_+^{(d-1)}(x)}{k!} - \\sum_{k=0}^{d-2} \\frac{(-x)^k \\phi^{(k)}(x)}{k!}
 ```
 """
-struct 𝒲₋₁{Tϕ,Tϕt} <: Distributions.ContinuousUnivariateDistribution
+function taylor(f, x, d, T)
+    return f(x + TaylorSeries.Taylor1(T,d)).coeffs
+end
+struct 𝒲₋₁{Tϕ} <: Distributions.ContinuousUnivariateDistribution
     ϕ::Tϕ
     d::Int64
     function 𝒲₋₁(ϕ,d)
-        return new{typeof(ϕ),typeof(ϕ_taylor)}(ϕ,d,ϕ_taylor)
+        @assert ϕ(0) == 1
+        @assert ϕ(Inf) == 0
+        # And assertion about d-monotony... how can this be check ? this is hard. 
+        return new{typeof(ϕ)}(ϕ,d)
     end
 end
 function Distributions.cdf(d::𝒲₋₁, x::Real)
-    rez = one(x)
-    t_taylor = TaylorSeries.Taylor1(Float64,d.d)
-    ϕ_taylor = d.ϕ(x + t_taylor).coeffs
-    ϕ_taylor[end] = max(ϕ_taylor[end], 0)
-    for k in 1:(d.d-1)
-        rez -= (-1)^k * x^k * ϕ_taylor[k+1]
+    rez = zero(x)
+    c_ϕ = taylor(d.ϕ, x, d.d, typeof(x))
+    c_ϕ[end] = max(c_ϕ[end], 0)
+    for k in 0:(d.d-1)
+        rez += (-1)^k * x^k * c_ϕ[k+1]
     end
-    return rez
+    return 1-rez
 end
 function Distributions.logpdf(d::𝒲₋₁, x::Real)
-    t_taylor = TaylorSeries.Taylor1(Float64,d.d+1)
-    ϕ_d = d.ϕ(x + t_taylor).coeffs[end]
-    return (d.d-1)*log(x) + log(ϕ_d) - sum(log.(1:(d.d-1)))
-
+    ϕ_d = taylor(d.ϕ, x, d.d+1, typeof(x))[end]
+    r = (d.d-1)*log(x) - sum(log.(1:(d.d-1)))
+    return log(ϕ_d) + r
 end
 function Distributions.rand(rng::Distributions.AbstractRNG, d::𝒲₋₁)
     u = rand(rng)

test/testing_the_paper.jl

@@ -1,136 +1,71 @@
-@testitem "Exemple 2.1" begin
+@testitem "Exemple 2.1 - WCopula, dimension 2" begin
     using Distributions
-    d=10
-    ϕ = 𝒲(Dirac(1),d)
-    generator_ex2_2(x) = max((1-x)^(d-1),0)
-    @test all(ϕ(x) ≈ generator_ex2_2(x) for x in 0:0.01:10)
+    d=2
+    X = Dirac(1)
+    ϕ(x, d) = max((1-x)^(d-1),zero(x))
+    Xhat = 𝒲₋₁(x -> ϕ(x,d),d)
+    ϕhat = 𝒲(X,d)
+
+    @test maximum(abs.([cdf(X,x) - cdf(Xhat,x) for x in 0:0.01:10*d])) <= sqrt(eps(Float64))
+    @test maximum(abs.([ϕ(x, d) - ϕhat(x) for x in 0:0.01:10])) <= sqrt(eps(Float64))
 end
 
-@testitem "Exemple 3.2" begin
+@testitem "Exemple 3.2 - IndependantCopula, dimension 10" begin
     using Distributions
-    d=10
-    ϕ = 𝒲(Erlang(10),d)
-    gen_indep(x) = exp(-x)
-    @test maximum(abs.([ϕ(x) - gen_indep(x) for x in 0:0.01:10])) <= sqrt(eps(Float64))
-    # @test all(ϕ(x) ≈ gen_indep(x) for x in 0:0.01:10)
+    for d in 3:20
+        X = Erlang(d)
+        ϕ(x) = exp(-x)
+        Xhat = 𝒲₋₁(ϕ,d)
+        ϕhat = 𝒲(X,d)
+
+        @test maximum(abs.([cdf(X,x) - cdf(Xhat,x) for x in 0:0.01:3*d])) <= sqrt(eps(Float64))
+        @test maximum(abs.([ϕ(x) - ϕhat(x) for x in 0:0.01:10])) <= sqrt(eps(Float64))
+    end
 end
 
-
-# @testitem "Exemple 3.3: inverse williamson clayton" begin
-#     using SpecialFunctions
+@testitem "Exemple 3.3: ClaytonCopula" begin
+    using SpecialFunctions, Distributions
 
-#     # exemple 3.3. : back to clayton. 
-#     gen_clayton(x,θ) = (1 + θ * x)^(-1/θ)
-#     function true_radial_cdf_for_clayton(x,θ,d)
-#         if x < 0
-#             return zero(x)
-#         end
-#         if θ < 0
-#             α = -1/θ
-#             if x >= α
-#                 return one(x)
-#             end
-#             rez = zero(x)
-#             θx = θ*x
-#             cst = log(-θx/(1+θx))
-#             @show x, cst
-#             for k in 0:(d-1)
-#                 rez += exp(loggamma(α+k+1) - loggamma(k+1) + k*cst)
-#             end
-#             rez *= (1+θx)^α/gamma(α+1)
-#             return 1-rez
-#         elseif θ == 0
-#             return exp(-x)
-#         else
-#             rez = zero(x)
-#             for k in 0:(d-1)
-#                 rez +=  prod(1+j*θ for j in 0:(k-1))/factorial(k) * x^k * (1+θ*x)^(-(1/θ+k))
-#             end
-#             return 1-rez
-#         end
-#     end
-#     θ = -0.3
-#     X = 𝒲₋₁(x -> gen_clayton(x,θ),2)
-
-#     @test maximum(abs.([true_radial_cdf_for_clayton(x,θ,2) - X.F(x) for x in 0:0.01:10])) <= sqrt(eps(Float64))
-# end
-
-
-
-
-
-
-# # An easy case:
-# F1(x) = (1 - exp(-x)) * (x > 0)
-# X = FromCDF(F1)
-# x = rand(X,10000)
-# Plots.plot(t -> StatsBase.ecdf(x)(t), 0, 10)
-# Plots.plot!(F1)
-
-# # A more involved one: 
-# F2(x) = 1*(x >= 2) # Dirac(1)
-# X = FromCDF(F2)
-# x = rand(X,10000)
-# Plots.plot(t -> StatsBase.ecdf(x)(t), 0, 4)
-# Plots.plot!(F2)
-
-# # A more involved one: 
-# F3(x) = Distributions.cdf(Distributions.Binomial(10,0.7),x)
-# X = FromCDF(F3)
-# x = rand(X,10000)
-# Plots.plot(t -> StatsBase.ecdf(x)(t), 0, 12)
-# Plots.plot!(F3)
-
-
-# # Final try: 
-# F4(x) = (F1(x)+F2(x)+F3(x) + x>=>)/3
-# X = FromCDF(F4)
-# x = rand(X,10000)
-# Plots.plot(t -> StatsBase.ecdf(x)(t), 0, 12)
-# Plots.plot!(F4)
-
-
-
-
-
-
-
-# # ϕ = 𝒲(Gamma(1,2),10)
-# ϕ = 𝒲(Dirac(1),10)
-# ϕ(0.0)
-# ϕ(Inf)
-# using Plots
-# plot(x -> ϕ(x),xlims=(0,10))
-
-
-# # Now we could implement the inverse case
-# # which allows to construct the random variable corresponding to a generator. 
-# # some precomputations might be necessary. 
-# # and a struct. 
-
-# # a certain generator from Ex 2.1: 
-# gen_ex22(x,d) = max((1-x)^(d-1),0)
-# # the clayton gen from ex 2.3: 
-# gen_clayton(x,d,θ) = (1 + θ * x)^(-1/θ)
-# # Note: θ = 0 generates the independence copula. 
-# # θ < 0 is interesting, as it is equal to gen_ex22 at the lower bound θ = -1/(d-1)
-
-# # example 3.1 : 
-# # X = Dirac(1) correspond to gen_ex22
-# ϕ = 𝒲(Dirac(1),10)
-# plot(x -> ϕ(x),xlims=(0,10))
-# plot!(x->gen_ex22(x,10))
-# # indeed same plot ! 
-
-
-# # exemple 3.2: 
-# # independent copula correspond to erlang distribution with parameter d
-# ϕ = 𝒲(Erlang(10),10)
-# gen_indep(x) = exp(-x)
-# plot(x -> ϕ(x),xlims=(0,10))
-# plot!(gen_indep)
-
-
-# we should check if a williamson d-transform of this distribution recovers the generator correctly. 
-
-
+    # exemple 3.3. : back to clayton. 
+    ϕ(x, θ) = max((1 + θ * x),zero(x))^(-1/θ)
+    function F(x, θ, d)
+        if x < 0
+            return zero(x)
+        end
+        α = -1/θ
+        if θ < 0
+            if x >= α
+                return one(x)
+            end
+            rez = zero(x)
+            x_α = x/α
+            for k in 0:(d-1)
+                rez += gamma(α+1)/gamma(α-k+1)/gamma(k+1) * (x_α)^k * (1 - x_α)^(α-k)
+            end
+            return 1-rez
+        elseif θ == 0
+            return exp(-x)
+        else
+            rez = zero(x)
+            for k in 0:(d-1)
+                pr = one(θ)
+                for j in 0:(k-1)
+                    pr *= (1+j*θ)
+                end
+                rez += pr / gamma(k+1) * x^k * (1 + θ * x)^(-(1/θ+k))
+            end
+            return 1-rez
+        end
+    end
+
+    for (d, θ) in (
+        (3, 1/7),
+        (2, -0.2), 
+        (10, -1/10),
+        (2, -1.0)
+    )
+        Xhat = 𝒲₋₁(x -> ϕ(x,θ),d)
+        @test maximum(abs.([F(x,θ,d) - cdf(Xhat,x) for x in 0:0.01:10])) <= sqrt(eps(Float64))
+    end
+
+end