diff --git a/docs/src/Liouville.md b/docs/src/Liouville.md
index 9b68c850..a8391593 100644
--- a/docs/src/Liouville.md
+++ b/docs/src/Liouville.md
@@ -4,17 +4,39 @@ CurrentModule = Copulas
 
 ## Liouville Copulas
 
-!!! note "Not merged yet !"
-    Liouville copulas are comming in this PR : https://github.com/lrnv/Copulas.jl/pull/83, but this is not merged yet. 
+Archimedean copulas have been widely used in the literature due to their nice decomposition properties and easy parametrization. The interested reader can refer to the extensive literature [hofert2010,hofert2013a,mcneil2010,cossette2017,cossette2018,genest2011a,dibernardino2013a,dibernardino2013a,dibernardino2016,cooray2018,spreeuw2014](@cite) on Archimedean copulas, their nesting extensions and most importantly their estimation. One major drawback of the Archimedean family is that these copulas have exchangeable marginals (i.e., $C(\bm u) = C(\mathrm{p}(\bm u))$ for any permutation $p(\bm u)$ of $u_1,...,u_d$): the dependence structure is symmetric, which might not be a wanted property. However, from the Radial-simplex expression, we can easily extrapolate a little and take for $\bm S$ a non-uniform distribution on the simplex. 
 
-Archimedean copulas have been widely used in the literature due to their nice decomposition properties and easy parametrization. The interested reader can refer to the extensive literature [hofert2010,hofert2013a,mcneil2010,cossette2017,cossette2018,genest2011a,dibernardino2013a,dibernardino2013a,dibernardino2016,cooray2018,spreeuw2014](@cite) on Archimedean copulas, their nesting extensions and most importantly their estimation. 
-
-One major drawback of the Archimedean family is that these copulas have exchangeable marginals (i.e., $C(\bm u) = C(\mathrm{p}(\bm u))$ for any permutation $p(\bm u)$ of $u_1,...,u_d$): the dependence structure is symmetric, which might not be a wanted property. However, from the Radial-simplex expression, we can easily extrapolate a little and take for $\bm S$ a non-uniform distribution on the simplex. 
+**Definition (Liouville Copulas):** For $R$ a positive random variable, $\bm\alpha \in N^d$ and $S \sim \mathrm{Dirichlet}(\bm \alpha)$ a Dirichlet random vector on the simplex, the copula of the random vector $R\bm S$ is called the Liouville copula with radial part $R$ and Simplex parameters $\alpha$. 
 
 Liouville's copulas share many properties with Archimedean copulas, but are not exchangeable anymore. This is an easy way to produce non-exchangeable dependence structures. See [cote2019](@cite) for a practical use of this property.
 
 Note that Dirichlet distributions are constructed as $\bm S = \frac{\bm G}{\langle \bm 1, \bm G\rangle}$, where $\bm G$ is a vector of independent Gamma distributions with unit scale (and potentially different shapes: taking all shapes equal yields the Archimedean case). 
 
+There are still a few properties that are interesting for the implementation, but wich requires a few notations. Let's denote by 
+
+$$\phi_{d}$$
+
+the inverse Williamson $d$-transform of $R$ for any integer $d$.
+
+Then the distribution function of the Liouville copula is given by 
+
+$$C(\bm u) = \sum_{\bm i \le \bm \alpha} \frac{(-1)^{\lvert\bm i\rvert}}{\bm i!} * \phi^{(\lvert\bm i\rvert)}(\lvert \bm x \rvert) * \bm x^{\bm i},$$
+
+where $x_i = F_{\phi_{\alpha_i}}^{-1}(u_i)$ for all i.
+
+And for sampling, the same kind of algorithm is availiable. 
+
+!!! note "Complexity of the matter"
+    Note that here we need access to the quantile functions of all the (inverse) Williamson $\alpha_i$-transform of $\phi_{\lvert\bm\alpha\rvert}$, which itself is the $\lvert\bm\alpha\rvert$-Williamson transfrom of $R$. This complexity adds up quickly, but the package internally uses fast Faa-di-bruno formula at compile time to overcome much of it. 
+
+You can sample and compute the cdf and pdf of *any* Liouville copula, even one for which you provide the generator and/or the Radial distribution yourself.
+
+
+```@docs
+LiouvilleCopula
+```
+
+
 ```@bibliography
 Pages = ["Liouville.md"]
 Canonical = false
diff --git a/src/Copulas.jl b/src/Copulas.jl
index a5348721..b71c2490 100644
--- a/src/Copulas.jl
+++ b/src/Copulas.jl
@@ -79,4 +79,7 @@ module Copulas
            MCopula,
            WCopula
 
+       include("LiouvilleCopula.jl")
+       export LiouvilleCopula
+
 end
diff --git a/src/LiouvilleCopula.jl b/src/LiouvilleCopula.jl
new file mode 100644
index 00000000..8ecdc167
--- /dev/null
+++ b/src/LiouvilleCopula.jl
@@ -0,0 +1,67 @@
+"""
+    LiouvilleCopula{d, T, TG}
+
+Fields: 
+* `α::NTuple{d,T}` : the weights for each dimension
+* `G::TG` : the generator <: Generator. 
+
+Constructor: 
+
+    LiouvilleCopula(α::Vector{T},G::Generator)
+
+The Liouville copula has a structure that resembles the [`ArchimedeanCopula`](@ref), when you look at it from it's radial-simplex decomposition. 
+
+Recalling that, for C an archimedean copula with generator ``\\phi``, if ``\\mathbf U \\sim C``, then ``U \\equal R \\mathbf S`` for a random vector ``\\mathbf S \\sim`` `Dirichlet(ones(d))`, that is uniformity on the d-variate simplex, and a non-negative random variable ``R`` that is the Williamson d-transform of `\\phi`. 
+
+The Liouville copula has exactly the same expression but using another Dirichlet distribution instead than uniformity over the simplex. 
+
+References: 
+* [mcneil2010](@cite) McNeil, A. J., & Nešlehová, J. (2010). From archimedean to liouville copulas. Journal of Multivariate Analysis, 101(8), 1772-1790.
+"""
+struct LiouvilleCopula{d,TG} <: Copula{d}
+    α::NTuple{d,Int}
+    G::TG
+    function LiouvilleCopula(α::Vector{Int},G::Generator)
+        d = length(α)
+        @assert sum(α) <= max_monotony(G) "The generator you provided is not monotonous enough (the monotony index must be greater than sum(α), and thus this copula does not exists."
+        return new{d, typeof(G)}(Tuple(α), G)
+    end
+end
+williamson_dist(C::LiouvilleCopula{d,TG}) where {d,TG} = williamson_dist(C.G,sum(C.α))
+
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T,CT<:LiouvilleCopula}
+    # By default, we use the williamson sampling. 
+    r = Distributions.rand(williamson_dist(C))
+    for i in 1:length(C)
+        x[i] = Distributions.rand(rng,Distributions.Gamma(C.α[i],1))
+        x[i] = x[i] * r
+        x[i] = Distributions.cdf(williamson_dist(C.G,C.α[i]),x[i])
+    end
+    return x
+end
+function _cdf(C::CT,v) where {CT<:LiouvilleCopula}
+    d = length(C)
+    if all(v == 1)
+        return one(eltype(v))
+    end
+    if any(v == 0)
+        return zero(eltype(v))
+    end
+    u = 1 .- v # Because survival.
+    Gis = [williamson_dist(C.G,C.α[i]) for i in 1:d]
+    x = [Distributions.quantile(Gis[i],u[i]) for i in 1:d]
+    sx = sum(x)
+    r = zero(eltype(u))
+    for i in CartesianIndices(C.α)
+        ii = (ij-1 for ij in Tuple(i))
+        sii = sum(ii)
+        fii = prod(factorial.(ii))
+        # This version is really not efficient as derivatives of the generator ϕ⁽ᵏ⁾(C.G, sii, sx) could be pre-computed since many sii will be the same and sx does never change. 
+        r += (-1)^sii / fii * ϕ⁽ᵏ⁾(C.G, sii, sx) * prod(x .^ii)
+    end
+
+    # r is the survival function, not the cdf ! 
+    # but if we use 1-u that might work. 
+
+    return r
+end
diff --git a/src/UnivariateDistribution/ClaytonWilliamsonDistribution.jl b/src/UnivariateDistribution/ClaytonWilliamsonDistribution.jl
index 6eaf3586..6de9e06a 100644
--- a/src/UnivariateDistribution/ClaytonWilliamsonDistribution.jl
+++ b/src/UnivariateDistribution/ClaytonWilliamsonDistribution.jl
@@ -33,7 +33,12 @@ function Distributions.cdf(D::ClaytonWilliamsonDistribution, x::Real)
         return 1-rez
     end
 end
-function Distributions.rand(rng::Distributions.AbstractRNG, d::ClaytonWilliamsonDistribution)
-    u = rand(rng)
+function _quantile(d::ClaytonWilliamsonDistribution,u)
     Roots.find_zero(x -> (Distributions.cdf(d,x) - u), (0.0, Inf))
+end
+function Distributions.rand(rng::Distributions.AbstractRNG, d::ClaytonWilliamsonDistribution)
+    _quantile(d,rand(rng))
+end
+function Distributions.quantile(d::ClaytonWilliamsonDistribution,p::Real)
+    _quantile(d,p)
 end
\ No newline at end of file
diff --git a/src/UnivariateDistribution/Logarithmic.jl b/src/UnivariateDistribution/Logarithmic.jl
index ac009355..24097534 100644
--- a/src/UnivariateDistribution/Logarithmic.jl
+++ b/src/UnivariateDistribution/Logarithmic.jl
@@ -1,6 +1,6 @@
 # Corresponds to https://en.wikipedia.org/wiki/Logarithmic_distribution
 struct Logarithmic{T<:Real} <: Distributions.DiscreteUnivariateDistribution
-    α::T # in (0,1), the weight of the logarithmic distribution.
+    α::T # in (0,1), the weight of the logarithmic distribution, 1-p of the wikipedia parameter p. 
     h::T # = ln(1-α), so in (-Inf,0)
     function Logarithmic(h::T) where {T <: Real}
         Tf = promote_type(T,Float64)
@@ -12,6 +12,31 @@ Base.eltype(::Logarithmic{T}) where T = promote_type(T,Float64)
 function Distributions.logpdf(d::Logarithmic{T}, x::Real) where T
     insupport(d, x) ? x*log1p(-d.α) - log(x) - log(-log(d.α)) : log(zero(T))
 end
+function _my_beta_inc(α,k)
+    # this computes \int_{0}^(1-α) t^k/(1-t) dt, for α in [0,1] and k a positive integer. 
+    r = zero(α)
+    pwα = 1
+    sgn = 1
+    for l in 0:k
+        r += sgn * binomial(k,l)  * (1 - pwα) / l
+        sgn *= -1
+        pwα *= α
+    end
+    return r
+end
+function Distributions.cdf(d::Logarithmic{T}, x::Real) where T
+    # Needs cdf to be implemented here. 
+    # quand even more quantile?
+    # ee her e: https://en.wikipedia.org/wiki/Logarithmic_distribution 
+    # return 1 + incbeta(p,k=1,0)/log1p(-p)
+    k = Int(trunc(x))
+    return 1 + _my_beta_inc(d.α,k)/log(d.α)
+    # SpecialFunctions.beta_inc(k+1,0,p)/log1p(-p) *SpecialFunctions.beta(k+1,0)
+end
+function Distributions.quantile(d::Logarithmic{T}, p::Real) where T
+    return Roots.find_zero(x -> Distributions.cdf(d,x) - p, (0, Inf))
+end
+
 function Distributions.rand(rng::Distributions.AbstractRNG, d::Logarithmic{T}) where T
     # Sample a Log(p) distribution with the algorithms "LK" and "LS" of Kemp (1981).
     if d.h > -3
diff --git a/src/UnivariateDistribution/WilliamsonFromFrailty.jl b/src/UnivariateDistribution/WilliamsonFromFrailty.jl
index 19288d7b..a18b0d7c 100644
--- a/src/UnivariateDistribution/WilliamsonFromFrailty.jl
+++ b/src/UnivariateDistribution/WilliamsonFromFrailty.jl
@@ -9,13 +9,27 @@ function Distributions.rand(rng::Distributions.AbstractRNG, D::WilliamsonFromFra
     sy = rand(rng,Distributions.Erlang(d))
     return sy/f
 end
-function Distributions.cdf(D::WilliamsonFromFrailty{TF,d}, x::Real) where {TF,d}
-    # how to compute this cdf ? 
-    return 1 - Distributions.expectation(
+function _p(D::WilliamsonFromFrailty{TF,d}, x) where {TF,d}
+
+    # Maybe we need to use a taylor serie approximation near x == 0
+    # since this is not working for x too small. 
+    # otherwise, we could simply truncate this: 
+    if x < sqrt(eps(eltype(x)))
+        return one(x)
+    end
+
+    return Distributions.expectation(
         e -> Distributions.cdf(D.frailty_dist,e/x),
         Distributions.Erlang(d)
     )
 end
+function Distributions.cdf(D::WilliamsonFromFrailty{TF,d}, x::Real) where {TF,d}
+    # how to compute this cdf ? 
+    return 1 - _p(D,x)
+end
+function Distributions.quantile(D::WilliamsonFromFrailty{TF,d},u::Real) where {TF,d}
+    return Roots.find_zero(x -> _p(D, x) - (1-u), (0, Inf))
+end
 function Distributions.pdf(D::WilliamsonFromFrailty{TF,d}, x::Real) where {TF,d}
     # how to compute this cdf ? 
     return 1/x^2 * Distributions.expectation(
diff --git a/test/liouville_tests.jl b/test/liouville_tests.jl
new file mode 100644
index 00000000..06fee840
--- /dev/null
+++ b/test/liouville_tests.jl
@@ -0,0 +1,24 @@
+
+@testitem "Test of Trivariate Liouville Copulas" begin
+    using Random, Distributions
+    using StableRNGs
+    using StatsBase
+    rng = StableRNG(123)
+
+    for G in (
+        Copulas.AMHGenerator(0.6),
+        Copulas.AMHGenerator(-0.3),
+        Copulas.ClaytonGenerator(-0.05),
+        Copulas.IndependentGenerator(),
+        Copulas.GumbelBarnettGenerator(0.7),
+        Copulas.InvGaussianGenerator(0.05),
+        Copulas.InvGaussianGenerator(8),
+        Copulas.WilliamsonGenerator(LogNormal(),20),
+    )
+        C = LiouvilleCopula(rand(1:6,3),G)
+        u = rand(C,10)
+        cdf(C,u)
+        # check margins uniformity ? 
+        @test true
+    end
+end
\ No newline at end of file
diff --git a/test/margins_uniformity.jl b/test/margins_uniformity.jl
index 48b6933f..945375f0 100644
--- a/test/margins_uniformity.jl
+++ b/test/margins_uniformity.jl
@@ -25,7 +25,18 @@
         FGMCopula(2,1),
         MCopula(4),
         PlackettCopula(2.0),
+        ArchimedeanCopula(3,WilliamsonGenerator(LogNormal(),4)),
         # Others ? Yes probably others too ! 
+        (LiouvilleCopula([1,5],G) for G in (
+            Copulas.AMHGenerator(0.6),
+            Copulas.AMHGenerator(-0.3),
+            Copulas.ClaytonGenerator(-0.05),
+            Copulas.IndependentGenerator(),
+            Copulas.GumbelBarnettGenerator(0.7),
+            Copulas.InvGaussianGenerator(0.05),
+            Copulas.InvGaussianGenerator(8),
+            Copulas.WilliamsonGenerator(LogNormal(),6),
+        ))...
     )
     n = 1000
     U = Uniform(0,1)