Merge pull request #136 from jojal5/dev

Dev

Project.toml

@@ -1,7 +1,7 @@
 name = "Extremes"
 uuid = "fe3fe864-1b39-11e9-20b8-1f96fa57382d"
 authors = ["Jonathan Jalbert <[email protected]>"]
-version = "1.0.2"
+version = "1.0.3"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

README.md

@@ -47,6 +47,20 @@ The datasets that are available through this package are the datasets referenced
 They were obtained using the R package `ismev`.  
 https://www.rdocumentation.org/packages/ismev/  
 
+## Lightweight fit functions
+
+The package also provides fit functions using maximum likelihood and probability-weighted moments for the Generalized Extreme Value, Gumbel, and Generalized Pareto distributions. They can be used with the `fit` function for the corresponding type in `Distributions.jl`. For example:
+
+```julia
+julia> using Distributions, Extremes
+
+y = rand(GeneralizedExtremeValue(0,1,.1))
+
+fit(GeneralizedExtremeValue, y)
+fit(GeneralizedExtremeValue, y, method="mle")
+fit(GeneralizedExtremeValue, y, method="pwm")
+```
+
 ## How to cite this library
 Jalbert, J., Farmer, M., Gobeil, G., & Roy, P. (2024). Extremes.jl: Extreme Value Analysis in Julia. Journal of Statistical Software, 109(6), 1–35. https://doi.org/10.18637/jss.v109.i06
 

ReleaseNotes.md

@@ -18,4 +18,7 @@
 - Refactor Hessian computations with PDMats.
 - Implement generic Delta method using PDMats.
 
+## 1.0.3
+- Add lightweight functions for simple fit of Generalized Extreme Value, Gumbel and Generalized Pareto distributions. See [`Extremes.fit`](@ref). 
+
 ## Nightly

src/Extremes.jl

@@ -9,7 +9,7 @@ using Gadfly
 import CSV
 import ForwardDiff, SpecialFunctions
 
-import Distributions: insupport, location, logpdf, maximum, minimum, params, pdf, quantile, scale, shape
+import Distributions: fit, fit_mle, insupport, location, logpdf, maximum, minimum, params, pdf, quantile, scale, shape
 import Statistics.var
 import Base.length, Base.maximum, Base.sum
 
@@ -44,6 +44,7 @@ export
     getcluster,
 
     # Fitting functions
+    fit,
     gevfit,
     gevfitbayes,
     gevfitpwm,

src/parameterestimation.jl

@@ -1,3 +1,5 @@
 include(joinpath("parameterestimation", "bayesian.jl"))
 include(joinpath("parameterestimation", "maximumlikelihood.jl"))
 include(joinpath("parameterestimation", "probabilityweightedmoment.jl"))
+include(joinpath("parameterestimation", "lightweight_functions.jl"))
+

src/parameterestimation/lightweight_functions.jl

@@ -0,0 +1,119 @@
+function fit_mle(pd::Type{<:GeneralizedExtremeValue}, y::Vector{<:Real})
+
+    fm = gevfit(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+function fit_pwm(pd::Type{<:GeneralizedExtremeValue}, y::Vector{<:Real})
+
+    fm = gevfitpwm(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+"""
+    fit(pd::Type{<:GeneralizedExtremeValue}, y::Vector{<:Real}; method::String="mle")
+
+Fit the Generalized Extreme Value distribution to the data vector `y`.
+
+## Details
+
+This function estimates the parameters of the GEV distribution from an identically distributed random sample. It is a lightweight feature of the Extremes.jl library.
+
+In the presence of non-stationarity or for advanced inference, the [`gevfit`](@ref) function is more appropriate.
+
+The `method` argument allows selecting the estimation method: `mle` for maximum likelihood estimation or `pwm` for the probability-weighted moments method.
+"""
+function fit(pd::Type{<:GeneralizedExtremeValue}, y::Vector{<:Real}; method::String="mle")
+    @assert method ∈ ("mle", "pwm") "Method $method is not a valid method among mle and pwm." 
+
+    if method == "mle"
+        return fit_mle(pd, y)
+    else
+        return fit_pwm(pd, y)
+    end
+
+end
+
+function fit_mle(pd::Type{<:Gumbel}, y::Vector{<:Real})
+
+    fm = gumbelfit(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+function fit_pwm(pd::Type{<:Gumbel}, y::Vector{<:Real})
+
+    fm = gumbelfitpwm(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+"""
+    fit(pd::Type{<:Gumbel}, y::Vector{<:Real}; method::String="mle")
+
+Fit the Gumbel distribution to the data vector `y`.
+
+## Details
+
+This function estimates the parameters of the Gumbel distribution from an identically distributed random sample. It is a lightweight feature of the Extremes.jl library.
+
+In the presence of non-stationarity or for advanced inference, the [`gevfit`](@ref) function is more appropriate.
+
+The `method` argument allows selecting the estimation method: `mle` for maximum likelihood estimation or `pwm` for the probability-weighted moments method.
+"""
+function fit(pd::Type{<:Gumbel}, y::Vector{<:Real}; method::String="mle")
+    @assert method ∈ ("mle", "pwm") "Method $method is not a valid method among mle and pwm." 
+
+    if method == "mle"
+        return fit_mle(pd, y)
+    else
+        return fit_pwm(pd, y)
+    end
+
+end
+
+function fit_mle(pd::Type{<:GeneralizedPareto}, y::Vector{<:Real})
+
+    fm = gpfit(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+function fit_pwm(pd::Type{<:GeneralizedPareto}, y::Vector{<:Real})
+
+    fm = gpfitpwm(y)
+    fd = Extremes.getdistribution(fm)[]
+
+    return fd
+end
+
+"""
+    fit(pd::Type{<:GeneralizedPareto}, y::Vector{<:Real}; method::String="mle")
+
+Fit the Generalized Pareto distribution to the exceedences vector `y`.
+
+## Details
+
+This function estimates the parameters of the GP distribution from an identically distributed random sample. It is a lightweight feature of the Extremes.jl library.
+
+In the presence of non-stationarity or for advanced inference, the [`gpfit`](@ref) function is more appropriate.
+
+The `method` argument allows selecting the estimation method: `mle` for maximum likelihood estimation or `pwm` for the probability-weighted moments method.
+"""
+function fit(pd::Type{<:GeneralizedPareto}, y::Vector{<:Real}; method::String="mle")
+    @assert method ∈ ("mle", "pwm") "Method $method is not a valid method among mle and pwm." 
+
+    if method == "mle"
+        return fit_mle(pd, y)
+    else
+        return fit_pwm(pd, y)
+    end
+
+end

test/parameterestimation/lightweight_functions_test.jl

@@ -0,0 +1,57 @@
+@testset "lightweight_functions.m" begin
+
+    @testset "fit(GeneralizedExtremeValue, y).jl" begin
+
+        df = CSV.read("dataset/gev_stationary.csv", DataFrame)
+        y = df.y
+
+        # Assert estimation method
+        @test_throws AssertionError fit(GeneralizedExtremeValue, y, method = "unsupported")
+
+        fm_mle = Extremes.fit_mle(GeneralizedExtremeValue, y)   
+
+        @test all(isapprox.(params(fm_mle),(0.0009, exp(0.0142), -0.0060), atol=.0001))
+
+        fm_pwm = Extremes.fit_pwm(GeneralizedExtremeValue, y) 
+
+        @test all(isapprox.(params(fm_pwm),(-0.0005, exp(0.0125), -0.0033), atol=.0001))
+
+    end
+
+    @testset "fit(Gumbel, y).jl" begin
+
+        df = CSV.read("dataset/gev_stationary.csv", DataFrame)
+        y = df.y
+
+        # Assert estimation method
+        @test_throws AssertionError fit(Gumbel, y, method = "unsupported")
+
+        fm_mle = Extremes.fit_mle(Gumbel, y)   
+
+        @test all(isapprox.(params(fm_mle),(-0.0023, exp(0.0124)), atol=.0001))
+
+        fm_pwm = Extremes.fit_pwm(Gumbel, y) 
+
+        @test all(isapprox.(params(fm_pwm),(-0.0020, exp(0.0095)), atol=.0001))
+
+    end
+
+    @testset "fit(GeneralizedPareto, y).jl" begin
+
+        df = CSV.read("dataset/gp_stationary.csv", DataFrame)
+        y = df.y
+
+        # Assert estimation method
+        @test_throws AssertionError fit(GeneralizedPareto, y, method = "unsupported")
+
+        fm_mle = Extremes.fit_mle(GeneralizedPareto, y)   
+
+        @test all(isapprox.(params(fm_mle), (0., exp(-0.0135), 0.0059), atol=.0001))
+
+        fm_pwm = Extremes.fit_pwm(GeneralizedPareto, y) 
+
+        @test all(isapprox.(params(fm_pwm),(0., exp(-0.0199), 0.0122), atol=.0001))
+
+    end
+
+end

test/parameterestimation_test.jl

@@ -2,5 +2,5 @@
     include(joinpath("parameterestimation", "bayesian_test.jl"))
     include(joinpath("parameterestimation", "maximumlikelihood_test.jl"))
     include(joinpath("parameterestimation", "probabilityweightedmoment_test.jl"))
-
+    include(joinpath("parameterestimation", "lightweight_functions_test.jl"))
 end

test/structures/AbstractFittedExtremeValueModel/bayesianAbstractExtremeValueModel_test.jl

@@ -24,7 +24,7 @@
     ci = cint(rl)
 
     @testset "getdistribution(fittedmodel)" begin
-        @test all(vec(Extremes.getdistribution(fm)) .== pd)
+        @test all(vec(Extremes.getdistribution(fm)) .≈ pd)
     end