fit Gamma by two quantiles

README.md

@@ -19,7 +19,7 @@ to
 
 This can also be used to approximate one distribution via a different distribution by matching its moments.
 
-User needs to [explicitly use Optim.jl](https://bgctw.github.io/DistributionFits.jl/stable/set_optimize/) for DitributionFits.jl to work properly:
+User needs to [explicitly using Optim.jl](https://bgctw.github.io/DistributionFits.jl/stable/set_optimize/) for DitributionFits.jl to work properly:
 ```julia
 use DistributionFits, Optim
 ```
@@ -29,5 +29,7 @@ Currently, support for the following distributios are implemented:
 - Lognormal
 - Logitnormal 
 - Exponential 
+- Weibull (only fitting to two quantiles)
+- Gamma (only fitting to two quantiles)
 
 See [Documentation](https://bgctw.github.io/DistributionFits.jl/dev)

docs/make.jl

@@ -21,6 +21,7 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Dependencies" => "set_optimize.md",
+        "Gamma" => "gamma.md",
         "LogNormal" => "lognormal.md",
         "LogitNormal" => "logitnormal.md",
         "Weibull" => "weibull.md",

docs/src/gamma.md

@@ -0,0 +1,18 @@
+```@meta
+CurrentModule = DistributionFits
+```
+
+# Gamma distribution
+
+The [Gamma distribution](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.Gamma) has a scale and a shape parameter.
+It can be fitted using two quantiles.
+
+```jldoctest; output = false, setup = :(using DistributionFits,Optim)
+d = fit(Gamma, @qp_m(0.5), @qp_uu(5));
+median(d) ≈ 0.5 && quantile(d, 0.975) ≈ 5
+# output
+true
+```
+
+Fitting by mean or mode is not yet implemented.
+

docs/src/weibull.md

@@ -4,7 +4,7 @@ CurrentModule = DistributionFits
 
 # Weibull distribution
 
-The Weibull distribution has a scale and a shape parameter.
+The [Weibull distribution](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.Weibull) has a scale and a shape parameter.
 It can be fitted using two quantiles.
 
 ```jldoctest; output = false, setup = :(using DistributionFits,Optim)

src/optimizer.jl

@@ -32,7 +32,7 @@ struct NotSetOptimizer <: AbstractDistributionFitOptimizer; end
 
 
 optimize(f, ::NotSetOptimizer, lower, upper) = error(
-    "Optimizer not set yet. Either invoke 'use Optim' or 'DistributionFits.set_optimizer(...)'.")
+    "Optimizer not set yet. Either invoke 'using Optim' or 'DistributionFits.set_optimizer(...)'.")
 
 
 optimizer = NotSetOptimizer();

src/univariate/continuous/gamma.jl

@@ -0,0 +1,43 @@
+# function fit(::Type{Gamma}, m::AbstractMoments)
+#     # https://en.wikipedia.org/wiki/Gamma_distribution
+#     n_moments(m) >= 1 || error("Need mean to estimate exponential")
+#     return Gamma(mean(m))
+# end
+
+function fit(::Type{Gamma}, lower::QuantilePoint, upper::QuantilePoint)
+    # https://www.johndcook.com/quantiles_parameters.pdf
+    #    Shape-scale paras, shape α - cook:α - w:k, scale θ - cook:β - w:θ
+    0 < lower.p < upper.p < 1 || error(
+        "Expected 0 < lower.p < upper.p < 1, but got " * 
+        "lower.p = $(lower.p) and upper.p = $(upper.p)")
+    0 < lower.q < upper.q || error(
+        "Expected 0 < lower.q < upper.q, but got " *
+        "lower.q = $(lower.q) and upper.q = $(upper.q)")
+    logp1 = log(lower.p)
+    logp2 = log(upper.p)
+    lhs = upper.q/lower.q
+    #α = 1
+    fcost = (α) -> begin
+        gamma1 = Gamma(α,1)
+        rhs = invlogcdf(gamma1, logp2)/invlogcdf(gamma1, logp1)
+        abs2(rhs - lhs)
+    end
+    resOpt = optimize(fcost, 1e-2, 1e8)
+    resOpt.converged || error("Could not fit distribution to quantiles. resOpt=$resOpt")
+    α_opt = resOpt.minimizer
+    fbeta(α) = lower.q / invlogcdf(Gamma(α,1), log(lower.p))
+    Gamma(α_opt, fbeta(α_opt))
+end
+
+# function fit_mean_quantile(::Type{Gamma}, mean::Real, qp::QuantilePoint)
+#     # only fit to mean
+#     fit(Type{Gamma}, AbstractMoments(mean))
+# end
+
+# function fit_mode_quantile(::Type{Gamma}, mode::Real, qp::QuantilePoint)
+#     # ignore mode (its always at 0)
+#     θ = -qp.q/log(1-qp.p)
+#     Gamma(θ)
+# end
+
+

src/univariates.jl

@@ -32,7 +32,8 @@ const continuous_distributions = [
      "exponential",
     # "fdist",
     # "frechet",
-    # "gamma", "erlang",
+    "gamma", 
+    # "erlang",
     # "pgeneralizedgaussian", # GeneralizedGaussian depends on Gamma
     # "generalizedpareto",
     # "generalizedextremevalue",

test/gamma.jl

@@ -0,0 +1,28 @@
+# @testset "fit moments" begin
+# end;
+
+# @testset "fit to quantilepoint and mode - ignore mode" begin
+# end;
+
+# @testset "fit two quantiles same" begin
+# end;
+
+@testset "fit two quantiles" begin
+    qpl = @qp_m(0.5)
+    qpu = @qp_uu(5)
+    d = fit(Gamma, qpl, qpu);
+    @test median(d) ≈ 0.5
+    @test quantile(d, 0.975) ≈ 5
+    #Plots.plot(d); 
+end;
+
+i_tmp = () -> begin
+    dw = fit(Weibull, @qp_m(0.5), @qp_uu(5));
+    dg = fit(Gamma, @qp_m(0.5), @qp_uu(5));
+    Plots.plot(dw)
+    Plots.plot!(dg, xlim=(0,0.2))
+    pdf(dw, 0)
+    pdf(dg, 0)
+    pdf(dw, 0.01)
+    pdf(dg, 0.01)
+end

test/runtests.jl

@@ -5,6 +5,10 @@ using Random: Random
 @testset "optimize error" begin
     @test_throws Exception DistributionFits.optimize(x -> x*x, DistributionFits.optimizer, -1, 1)
 end
+# Optim package for interactive testing
+i_loadlibs = () -> begin
+    push!(LOAD_PATH, expanduser("~/julia/scimltools/")) # access local package repo
+end
 using Optim: Optim, optimize
 @testset "optimize set in __init__ after using Optim" begin
     # set in __init__
@@ -18,6 +22,7 @@ const tests = [
     "logitnormal",
     "exponential",
     "weibull",
+    "gamma",
 ]
 #tests = ["logitnormal"]