/
pareto.jl
132 lines (97 loc) · 3.59 KB
/
pareto.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
"""
Pareto(α,θ)
The *Pareto distribution* with shape `α` and scale `θ` has probability density function
```math
f(x; \\alpha, \\theta) = \\frac{\\alpha \\theta^\\alpha}{x^{\\alpha + 1}}, \\quad x \\ge \\theta
```
```julia
Pareto() # Pareto distribution with unit shape and unit scale, i.e. Pareto(1, 1)
Pareto(α) # Pareto distribution with shape α and unit scale, i.e. Pareto(α, 1)
Pareto(α, θ) # Pareto distribution with shape α and scale θ
params(d) # Get the parameters, i.e. (α, θ)
shape(d) # Get the shape parameter, i.e. α
scale(d) # Get the scale parameter, i.e. θ
```
External links
* [Pareto distribution on Wikipedia](http://en.wikipedia.org/wiki/Pareto_distribution)
"""
struct Pareto{T<:Real} <: ContinuousUnivariateDistribution
α::T
θ::T
Pareto{T}(α::T, θ::T) where {T} = new{T}(α, θ)
end
function Pareto(α::T, θ::T; check_args::Bool=true) where {T <: Real}
@check_args Pareto (α, α > zero(α)) (θ, θ > zero(θ))
return Pareto{T}(α, θ)
end
Pareto(α::Real, θ::Real; check_args::Bool=true) = Pareto(promote(α, θ)...; check_args=check_args)
Pareto(α::Integer, θ::Integer; check_args::Bool=true) = Pareto(float(α), float(θ); check_args=check_args)
Pareto(α::Real; check_args::Bool=true) = Pareto(α, one(α); check_args=check_args)
Pareto() = Pareto{Float64}(1.0, 1.0)
@distr_support Pareto d.θ Inf
#### Conversions
convert(::Type{Pareto{T}}, α::Real, θ::Real) where {T<:Real} = Pareto(T(α), T(θ))
Base.convert(::Type{Pareto{T}}, d::Pareto) where {T<:Real} = Pareto{T}(T(d.α), T(d.θ))
Base.convert(::Type{Pareto{T}}, d::Pareto{T}) where {T<:Real} = d
#### Parameters
shape(d::Pareto) = d.α
scale(d::Pareto) = d.θ
params(d::Pareto) = (d.α, d.θ)
@inline partype(d::Pareto{T}) where {T<:Real} = T
#### Statistics
function mean(d::Pareto{T}) where T<:Real
(α, θ) = params(d)
α > 1 ? α * θ / (α - 1) : T(Inf)
end
median(d::Pareto) = ((α, θ) = params(d); θ * 2^(1/α))
mode(d::Pareto) = d.θ
function var(d::Pareto{T}) where T<:Real
(α, θ) = params(d)
α > 2 ? (θ^2 * α) / ((α - 1)^2 * (α - 2)) : T(Inf)
end
function skewness(d::Pareto{T}) where T<:Real
α = shape(d)
α > 3 ? ((2(1 + α)) / (α - 3)) * sqrt((α - 2) / α) : T(NaN)
end
function kurtosis(d::Pareto{T}) where T<:Real
α = shape(d)
α > 4 ? (6(α^3 + α^2 - 6α - 2)) / (α * (α - 3) * (α - 4)) : T(NaN)
end
entropy(d::Pareto) = ((α, θ) = params(d); log(θ / α) + 1 / α + 1)
#### Evaluation
function pdf(d::Pareto{T}, x::Real) where T<:Real
(α, θ) = params(d)
x >= θ ? α * (θ / x)^α * (1/x) : zero(T)
end
function logpdf(d::Pareto{T}, x::Real) where T<:Real
(α, θ) = params(d)
x >= θ ? log(α) + α * log(θ) - (α + 1) * log(x) : -T(Inf)
end
function ccdf(d::Pareto, x::Real)
α, θ = params(d)
return (θ / max(x, θ))^α
end
cdf(d::Pareto, x::Real) = 1 - ccdf(d, x)
function logccdf(d::Pareto, x::Real)
α, θ = params(d)
return xlogy(α, θ / max(x, θ))
end
logcdf(d::Pareto, x::Real) = log1mexp(logccdf(d, x))
cquantile(d::Pareto, p::Real) = d.θ / p^(1 / d.α)
quantile(d::Pareto, p::Real) = cquantile(d, 1 - p)
#### Sampling
rand(rng::AbstractRNG, d::Pareto) = d.θ * exp(randexp(rng) / d.α)
## Fitting
function fit_mle(::Type{<:Pareto}, x::AbstractArray{T}) where T<:Real
# Based on
# https://en.wikipedia.org/wiki/Pareto_distribution#Parameter_estimation
θ = minimum(x)
n = length(x)
lθ = log(θ)
temp1 = zero(T)
for i=1:n
temp1 += log(x[i]) - lθ
end
α = n/temp1
return Pareto(α, θ)
end