Implement fit_mle for pareto distribution

Ph0non · Aug 20, 2015 · 1b4f4cb · 1b4f4cb
1 parent 62df2e8
commit 1b4f4cb
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 0 deletions.
diff --git a/doc/source/fit.rst b/doc/source/fit.rst
@@ -49,6 +49,7 @@ The ``fit_mle`` method has been implemented for the following distributions:
 - :ref:`gamma`
 - :ref:`geometric`
 - :ref:`laplace`
+- :ref:`pareto`
 - :ref:`poisson`
 - :ref:`uniform`
 

diff --git a/src/univariate/continuous/pareto.jl b/src/univariate/continuous/pareto.jl
@@ -78,3 +78,23 @@ quantile(d::Pareto, p::Float64) = cquantile(d, 1.0 - p)
 #### Sampling
 
 rand(d::Pareto) = d.θ * exp(randexp() / d.α)
+
+
+## Fitting
+
+function fit_mle{T <: Real}(::Type{Pareto}, x::AbstractArray{T})
+    # 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
diff --git a/test/fit.jl b/test/fit.jl
@@ -318,6 +318,14 @@ d = fit(Laplace, rand(Laplace(5.0, 3.0), N))
 @test_approx_eq_eps location(d) 5.0 0.1
 @test_approx_eq_eps scale(d) 3.0 0.2
 
+# Pareto
+
+x = rand(Pareto(3., 7.), n0)
+d = fit(Pareto, x)
+
+@test isa(d, Pareto)
+@test_approx_eq_eps shape(d) 3. 0.1
+@test_approx_eq_eps scale(d) 7. 0.1
 
 # Poisson