Skip to content

Commit bbc8a8b

Browse files
ararslanararslan
committed
Add the Kumaraswamy distribution
1 parent ec68da3 commit bbc8a8b

File tree

9 files changed

+384
-1
lines changed

9 files changed

+384
-1
lines changed

‎docs/src/univariate.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,13 @@ KSDist
284284
KSOneSided
285285
```
286286

287+
```@docs
288+
Kumaraswamy
289+
```
290+
```@example plotdensity
291+
plotdensity((0, 1), Kumaraswamy, (2, 5)) # hide
292+
```
293+
287294
```@docs
288295
Laplace
289296
```

‎src/Distributions.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ import StatsBase: kurtosis, skewness, entropy, mode, modes,
2424
import PDMats: dim, PDMat, invquad
2525

2626
using SpecialFunctions
27+
using Base.MathConstants
2728

2829
export
2930
# re-export Statistics
@@ -108,6 +109,7 @@ export
108109
Kolmogorov,
109110
KSDist,
110111
KSOneSided,
112+
Kumaraswamy,
111113
Laplace,
112114
Levy,
113115
Lindley,
@@ -345,7 +347,8 @@ Supported distributions:
345347
Frechet, FullNormal, FullNormalCanon, Gamma, GeneralizedPareto,
346348
GeneralizedExtremeValue, Geometric, Gumbel, Hypergeometric,
347349
InverseWishart, InverseGamma, InverseGaussian, IsoNormal,
348-
IsoNormalCanon, Kolmogorov, KSDist, KSOneSided, Laplace, Levy, Lindley, LKJ, LKJCholesky,
350+
IsoNormalCanon, Kolmogorov, KSDist, KSOneSided, Kumaraswamy,
351+
Laplace, Levy, Lindley, LKJ, LKJCholesky,
349352
Logistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal,
350353
MatrixTDist, MixtureModel, Multinomial,
351354
MultivariateNormal, MvLogNormal, MvNormal, MvNormalCanon,

‎src/univariate/continuous/kumaraswamy.jl

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
"""
2+
Kumaraswamy(a, b)
3+
4+
The *Kumaraswamy distribution* with shape parameters `a > 0` and `b > 0` has probability
5+
density function
6+
7+
```math
8+
f(x; a, b) = a b x^{a - 1} (1 - x^a)^{b - 1}, \\quad 0 < x < 1
9+
```
10+
11+
It is related to the [Beta distribution](@ref Beta) by the following identity:
12+
if ``X \\sim \\operatorname{Kumaraswamy}(a, b)`` then ``X^a \\sim \\operatorname{Beta}(1, b)``.
13+
Additionally, if ``X \\sim \\operatorname{Kumaraswamy}(1, 1)`` then
14+
``X \\sim \\operatorname{Uniform}(0, 1)``.
15+
16+
External links
17+
18+
- [Kumaraswamy distribution on Wikipedia](https://en.wikipedia.org/wiki/Kumaraswamy_distribution)
19+
20+
References
21+
22+
- Kumaraswamy, P. (1980). A generalized probability density function for double-bounded
23+
random processes. Journal of Hydrology. 46(1-2), 79-88.
24+
"""
25+
struct Kumaraswamy{T<:Real} <: ContinuousUnivariateDistribution
26+
a::T
27+
b::T
28+
end
29+
30+
function Kumaraswamy(a::Real, b::Real; check_args::Bool=true)
31+
@check_args Kumaraswamy (a, a > zero(a)) (b, b > zero(b))
32+
a′, b′ = promote(a, b)
33+
return Kumaraswamy{typeof(a′)}(a′, b′)
34+
end
35+
36+
Kumaraswamy() = Kumaraswamy{Float64}(1.0, 1.0)
37+
38+
Base.convert(::Type{Kumaraswamy{T}}, d::Kumaraswamy) where {T} = Kumaraswamy{T}(T(d.a), T(d.b))
39+
Base.convert(::Type{Kumaraswamy{T}}, d::Kumaraswamy{T}) where {T} = d
40+
41+
@distr_support Kumaraswamy 0 1
42+
43+
### Parameters
44+
45+
params(d::Kumaraswamy) = (d.a, d.b)
46+
partype(::Kumaraswamy{T}) where {T} = T
47+
48+
### Evaluation
49+
50+
function pdf(d::Kumaraswamy, x::Real)
51+
a, b = params(d)
52+
y = a * b * x^(a - 1) * (1 - x^a)^(b - 1)
53+
return 0 < x < 1 ? y : (isnan(x) ? oftype(y, NaN) : zero(y))
54+
end
55+
56+
function logpdf(d::Kumaraswamy, x::Real)
57+
a, b = params(d)
58+
_x = clamp(x, 0, 1) # Ensures we can still get a value when outside the support
59+
y = log(a) + log(b) + (a - 1) * log(_x) + (b - 1) * log1p(-_x^a)
60+
return 0 < x < 1 ? y : (isnan(x) ? oftype(y, NaN) : oftype(y, -Inf))
61+
end
62+
63+
function ccdf(d::Kumaraswamy, x::Real)
64+
a, b = params(d)
65+
y = (1 - clamp(x, 0, 1)^a)^b
66+
return x < 0 ? one(y) : (x > 1 ? zero(y) : y)
67+
end
68+
69+
cdf(d::Kumaraswamy, x::Real) = 1 - ccdf(d, x)
70+
71+
function logccdf(d::Kumaraswamy, x::Real)
72+
a, b = params(d)
73+
y = b * log1p(-clamp(x, 0, 1)^a)
74+
return x < 0 ? zero(y) : (x > 1 ? oftype(y, -Inf) : y)
75+
end
76+
77+
logcdf(d::Kumaraswamy, x::Real) = log1mexp(logccdf(d, x))
78+
79+
function quantile(d::Kumaraswamy, q::Real)
80+
a, b = params(d)
81+
return (1 - (1 - q)^inv(b))^inv(a)
82+
end
83+
84+
function entropy(d::Kumaraswamy)
85+
a, b = params(d)
86+
H = digamma(b + 1) + eulergamma
87+
return (1 - inv(b)) + (1 - inv(a)) * H - log(a) - log(b)
88+
end
89+
90+
function gradlogpdf(d::Kumaraswamy, x::Real)
91+
a, b = params(d)
92+
_x = clamp(x, 0, 1)
93+
_xᵃ = _x^a
94+
y = (a * (b * _xᵃ - 1) + (1 - _xᵃ)) / (_x * (_xᵃ - 1))
95+
return 0 < x < 1 ? y : (isnan(x) ? oftype(y, NaN) : zero(y))
96+
end
97+
98+
### Sampling
99+
100+
rand(rng::AbstractRNG, d::Kumaraswamy) = quantile(d, rand(rng))
101+
102+
### Statistics
103+
104+
_kumomentaswamy(a, b, n) = b * beta(1 + n / a, b)
105+
106+
mean(d::Kumaraswamy) = _kumomentaswamy(params(d)..., 1)
107+
108+
function var(d::Kumaraswamy)
109+
a, b = params(d)
110+
m₁ = _kumomentaswamy(a, b, 1)
111+
m₂ = _kumomentaswamy(a, b, 2)
112+
return m₂ - m₁^2
113+
end
114+
115+
function skewness(d::Kumaraswamy)
116+
a, b = params(d)
117+
μ = mean(d)
118+
σ² = var(d)
119+
m₂ = _kumomentaswamy(a, b, 2)
120+
m₃ = _kumomentaswamy(a, b, 3)
121+
return (2m₃ - μ * (3m₂ - μ^2)) / (σ² * sqrt(σ²))
122+
end
123+
124+
function kurtosis(d::Kumaraswamy)
125+
a, b = params(d)
126+
μ = mean(d)
127+
m₂ = _kumomentaswamy(a, b, 2)
128+
m₃ = _kumomentaswamy(a, b, 3)
129+
m₄ = _kumomentaswamy(a, b, 4)
130+
return (m₄ + μ * (-4m₃ + μ * (6m₂ - 3μ^2))) / var(d)^2 - 3
131+
end
132+
133+
median(d::Kumaraswamy) = (1 - 2^-inv(b))^inv(a)
134+
135+
function mode(d::Kumaraswamy)
136+
a, b = params(d)
137+
m = ((a - 1) / (a * b - 1))^inv(a)
138+
return a >= 1 && b >= 1 && !(a == b == 1) ? m : oftype(m, NaN)
139+
end

‎src/univariates.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -690,6 +690,7 @@ const continuous_distributions = [
690690
"kolmogorov",
691691
"ksdist",
692692
"ksonesided",
693+
"kumaraswamy",
693694
"laplace",
694695
"levy",
695696
"lindley",

‎test/ref/continuous/kumaraswamy.R

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
Kumaraswamy <- R6Class("Kumaraswamy",
2+
inherit=ContinuousDistribution,
3+
public=list(names=c("a", "b"),
4+
a=NA,
5+
b=NA,
6+
initialize=function(a=1, b=1) {
7+
self$a <- a
8+
self$b <- b
9+
},
10+
supp=function() { c(0, 1) },
11+
properties=function() { list() },
12+
pdf=function(x, log=FALSE) { dkumar(x, self$a, self$b, log=log) },
13+
cdf=function(x) { pkumar(x, self$a, self$b) },
14+
quan=function(x) { qkumar(x, self$a, self$b) }))

‎test/ref/continuous_test.lst

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,13 @@ InverseGaussian(1.0, 1.0)
9090
InverseGaussian(2.0, 1.5)
9191
InverseGaussian(2.0, 7.0)
9292

93+
Kumaraswamy()
94+
Kumaraswamy(0.5, 0.5)
95+
Kumaraswamy(5, 1.0)
96+
Kumaraswamy(1.0, 3)
97+
Kumaraswamy(2, 2)
98+
Kumaraswamy(2, 5)
99+
93100
Laplace()
94101
Laplace(2.0)
95102
Laplace(0.0, 1.0)

‎test/ref/continuous_test.ref.json

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2530,6 +2530,162 @@
25302530
{ "q": 0.90, "x": 3.38940701260244 }
25312531
]
25322532
},
2533+
{
2534+
"expr": "Kumaraswamy()",
2535+
"dtype": "Kumaraswamy",
2536+
"minimum": 0,
2537+
"maximum": 1,
2538+
"properties": {
2539+
},
2540+
"points": [
2541+
{ "x": 0.1, "pdf": 1, "logpdf": 0, "cdf": 0.1 },
2542+
{ "x": 0.2, "pdf": 1, "logpdf": 0, "cdf": 0.2 },
2543+
{ "x": 0.3, "pdf": 1, "logpdf": 0, "cdf": 0.3 },
2544+
{ "x": 0.4, "pdf": 1, "logpdf": 0, "cdf": 0.4 },
2545+
{ "x": 0.5, "pdf": 1, "logpdf": 0, "cdf": 0.5 },
2546+
{ "x": 0.6, "pdf": 1, "logpdf": 0, "cdf": 0.6 },
2547+
{ "x": 0.7, "pdf": 1, "logpdf": 0, "cdf": 0.7 },
2548+
{ "x": 0.8, "pdf": 1, "logpdf": 0, "cdf": 0.8 },
2549+
{ "x": 0.9, "pdf": 1, "logpdf": 0, "cdf": 0.9 }
2550+
],
2551+
"quans": [
2552+
{ "q": 0.10, "x": 0.1 },
2553+
{ "q": 0.25, "x": 0.25 },
2554+
{ "q": 0.50, "x": 0.5 },
2555+
{ "q": 0.75, "x": 0.75 },
2556+
{ "q": 0.90, "x": 0.9 }
2557+
]
2558+
},
2559+
{
2560+
"expr": "Kumaraswamy(0.5, 0.5)",
2561+
"dtype": "Kumaraswamy",
2562+
"minimum": 0,
2563+
"maximum": 1,
2564+
"properties": {
2565+
},
2566+
"points": [
2567+
{ "x": 0.0361, "pdf": 1.46198830409357, "logpdf": 0.379797361359587, "cdf": 0.1 },
2568+
{ "x": 0.1296, "pdf": 0.868055555555556, "logpdf": -0.141499562273699, "cdf": 0.2 },
2569+
{ "x": 0.2601, "pdf": 0.700280112044818, "logpdf": -0.356274863917393, "cdf": 0.3 },
2570+
{ "x": 0.4096, "pdf": 0.651041666666667, "logpdf": -0.42918163472548, "cdf": 0.4 },
2571+
{ "x": 0.5625, "pdf": 0.666666666666667, "logpdf": -0.405465108108164, "cdf": 0.5 },
2572+
{ "x": 0.7056, "pdf": 0.744047619047619, "logpdf": -0.295650242100958, "cdf": 0.6 },
2573+
{ "x": 0.8281, "pdf": 0.915750915750916, "logpdf": -0.0880108773227132, "cdf": 0.7 },
2574+
{ "x": 0.9216, "pdf": 1.30208333333333, "logpdf": 0.263965545834465, "cdf": 0.8 },
2575+
{ "x": 0.9801, "pdf": 2.52525252525252, "logpdf": 0.926341067727656, "cdf": 0.9 }
2576+
],
2577+
"quans": [
2578+
{ "q": 0.10, "x": 0.0361 },
2579+
{ "q": 0.25, "x": 0.19140625 },
2580+
{ "q": 0.50, "x": 0.5625 },
2581+
{ "q": 0.75, "x": 0.87890625 },
2582+
{ "q": 0.90, "x": 0.9801 }
2583+
]
2584+
},
2585+
{
2586+
"expr": "Kumaraswamy(5, 1.0)",
2587+
"dtype": "Kumaraswamy",
2588+
"minimum": 0,
2589+
"maximum": 1,
2590+
"properties": {
2591+
},
2592+
"points": [
2593+
{ "x": 0.630957344480193, "pdf": 0.792446596230557, "logpdf": -0.232630161961136, "cdf": 0.1 },
2594+
{ "x": 0.724779663677696, "pdf": 1.37972966146121, "logpdf": 0.32188758248682, "cdf": 0.2 },
2595+
{ "x": 0.786003085596623, "pdf": 1.90838945480909, "logpdf": 0.646259668973352, "cdf": 0.3 },
2596+
{ "x": 0.832553207401873, "pdf": 2.40224886796286, "logpdf": 0.876405326934776, "cdf": 0.4 },
2597+
{ "x": 0.870550563296124, "pdf": 2.87174588749259, "logpdf": 1.05492016798614, "cdf": 0.5 },
2598+
{ "x": 0.902880451447434, "pdf": 3.32269902974487, "logpdf": 1.20077741342131, "cdf": 0.6 },
2599+
{ "x": 0.931149915094838, "pdf": 3.75879323325023, "logpdf": 1.32409795728311, "cdf": 0.7 },
2600+
{ "x": 0.956352499790037, "pdf": 4.18255821036509, "logpdf": 1.43092307138273, "cdf": 0.8 },
2601+
{ "x": 0.979148362360977, "pdf": 4.59583059420061, "logpdf": 1.52514949990784, "cdf": 0.9 }
2602+
],
2603+
"quans": [
2604+
{ "q": 0.10, "x": 0.630957344480193 },
2605+
{ "q": 0.25, "x": 0.757858283255199 },
2606+
{ "q": 0.50, "x": 0.870550563296124 },
2607+
{ "q": 0.75, "x": 0.944087511294902 },
2608+
{ "q": 0.90, "x": 0.979148362360977 }
2609+
]
2610+
},
2611+
{
2612+
"expr": "Kumaraswamy(1.0, 3)",
2613+
"dtype": "Kumaraswamy",
2614+
"minimum": 0,
2615+
"maximum": 1,
2616+
"properties": {
2617+
},
2618+
"points": [
2619+
{ "x": 0.0345106153943703, "pdf": 2.79650925535847, "logpdf": 1.02837194489623, "cdf": 0.1 },
2620+
{ "x": 0.0716822332774442, "pdf": 2.58532162803826, "logpdf": 0.949849921125303, "cdf": 0.2 },
2621+
{ "x": 0.112095998257399, "pdf": 2.36512054893157, "logpdf": 0.860828992708955, "cdf": 0.3 },
2622+
{ "x": 0.156567334698251, "pdf": 2.13413598269404, "logpdf": 0.758061872824116, "cdf": 0.4 },
2623+
{ "x": 0.2062994740159, "pdf": 1.88988157484231, "logpdf": 0.636514168294813, "cdf": 0.5 },
2624+
{ "x": 0.263193700271923, "pdf": 1.62865056995694, "logpdf": 0.487751800752006, "cdf": 0.6 },
2625+
{ "x": 0.33056704991783, "pdf": 1.34442142396715, "logpdf": 0.295963752450819, "cdf": 0.7 },
2626+
{ "x": 0.415196452357427, "pdf": 1.02598556800602, "logpdf": 0.0256536803787092, "cdf": 0.8 },
2627+
{ "x": 0.535841116638722, "pdf": 0.646330407009565, "logpdf": -0.436444439994588, "cdf": 0.9 }
2628+
],
2629+
"quans": [
2630+
{ "q": 0.10, "x": 0.0345106153943703 },
2631+
{ "q": 0.25, "x": 0.0914397035839302 },
2632+
{ "q": 0.50, "x": 0.2062994740159 },
2633+
{ "q": 0.75, "x": 0.370039475052563 },
2634+
{ "q": 0.90, "x": 0.535841116638722 }
2635+
]
2636+
},
2637+
{
2638+
"expr": "Kumaraswamy(2, 2)",
2639+
"dtype": "Kumaraswamy",
2640+
"minimum": 0,
2641+
"maximum": 1,
2642+
"properties": {
2643+
},
2644+
"points": [
2645+
{ "x": 0.226531900511796, "pdf": 0.859628121964726, "logpdf": -0.151255399573567, "cdf": 0.1 },
2646+
{ "x": 0.324919696232906, "pdf": 1.16246804480858, "logpdf": 0.150545369764855, "cdf": 0.2 },
2647+
{ "x": 0.40415340338283, "pdf": 1.35255598879246, "logpdf": 0.301996127401067, "cdf": 0.3 },
2648+
{ "x": 0.47476660661689, "pdf": 1.4710105286101, "logpdf": 0.385949599044466, "cdf": 0.4 },
2649+
{ "x": 0.541196100146197, "pdf": 1.53073372946036, "logpdf": 0.42574718219016, "cdf": 0.5 },
2650+
{ "x": 0.606254458100165, "pdf": 1.53371594338211, "logpdf": 0.427693511992613, "cdf": 0.6 },
2651+
{ "x": 0.672515756317154, "pdf": 1.47340820005021, "logpdf": 0.387578220644201, "cdf": 0.7 },
2652+
{ "x": 0.743496068920369, "pdf": 1.33000620088785, "logpdf": 0.285183604544486, "cdf": 0.8 },
2653+
{ "x": 0.82690521463053, "pdf": 1.04596155492114, "logpdf": 0.0449366105897816, "cdf": 0.9 }
2654+
],
2655+
"quans": [
2656+
{ "q": 0.10, "x": 0.226531900511796 },
2657+
{ "q": 0.25, "x": 0.366025403784439 },
2658+
{ "q": 0.50, "x": 0.541196100146197 },
2659+
{ "q": 0.75, "x": 0.707106781186548 },
2660+
{ "q": 0.90, "x": 0.82690521463053 }
2661+
]
2662+
},
2663+
{
2664+
"expr": "Kumaraswamy(2, 5)",
2665+
"dtype": "Kumaraswamy",
2666+
"minimum": 0,
2667+
"maximum": 1,
2668+
"properties": {
2669+
},
2670+
"points": [
2671+
{ "x": 0.144400961350758, "pdf": 1.32728471201559, "logpdf": 0.28313528547594, "cdf": 0.1 },
2672+
{ "x": 0.208919841589934, "pdf": 1.7476387975003, "logpdf": 0.558265618295218, "cdf": 0.2 },
2673+
{ "x": 0.262392997058158, "pdf": 1.9725620435889, "logpdf": 0.679333227534167, "cdf": 0.3 },
2674+
{ "x": 0.311640094584387, "pdf": 2.07097247981028, "logpdf": 0.728018293967013, "cdf": 0.4 },
2675+
{ "x": 0.359790823540395, "pdf": 2.0664556357194, "logpdf": 0.725834886420501, "cdf": 0.5 },
2676+
{ "x": 0.409202630243412, "pdf": 1.96601311053933, "logpdf": 0.676007690339223, "cdf": 0.6 },
2677+
{ "x": 0.462598005187417, "pdf": 1.76563430983077, "logpdf": 0.568510008163275, "cdf": 0.7 },
2678+
{ "x": 0.524614464461574, "pdf": 1.44765227489845, "logpdf": 0.369943123480185, "cdf": 0.8 },
2679+
{ "x": 0.607488811024373, "pdf": 0.962804881088825, "logpdf": -0.0379045034051143, "cdf": 0.9 }
2680+
],
2681+
"quans": [
2682+
{ "q": 0.10, "x": 0.144400961350758 },
2683+
{ "q": 0.25, "x": 0.236458217673013 },
2684+
{ "q": 0.50, "x": 0.359790823540395 },
2685+
{ "q": 0.75, "x": 0.492078974093388 },
2686+
{ "q": 0.90, "x": 0.607488811024373 }
2687+
]
2688+
},
25332689
{
25342690
"expr": "Laplace()",
25352691
"dtype": "Laplace",

‎test/ref/rdistributions.R

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ source("continuous/generalizedpareto.R")
5656
source("continuous/gumbel.R")
5757
source("continuous/inversegamma.R")
5858
source("continuous/inversegaussian.R")
59+
source("continuous/kumaraswamy.R")
5960
source("continuous/laplace.R")
6061
source("continuous/levy.R")
6162
source("continuous/lindley.R")

0 commit comments

Comments
 (0)