From b4013855c9bab8604b9533bce83b768cfd6936d7 Mon Sep 17 00:00:00 2001
From: Oskar Laverny <oskar.laverny@gmail.com>
Date: Fri, 24 Nov 2023 15:57:02 +0100
Subject: [PATCH] Add functionality to the logarithmic distributons

---
 src/UnivariateDistribution/Logarithmic.jl | 27 ++++++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

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