From 0dac0778ebbfadcbdccc20f614ecc3e0c663b184 Mon Sep 17 00:00:00 2001
From: Thomas Wutzler <twutz@bgc-jena.mpg.de>
Date: Mon, 16 May 2022 13:29:14 +0200
Subject: [PATCH] fit Weibull by two quantiles

---
 docs/make.jl                             |  1 +
 docs/src/weibull.md                      | 18 ++++++++++++++++
 src/univariate/continuous/exponential.jl |  4 ++--
 src/univariate/continuous/weibull.jl     | 25 ++++++++++++++++++++++
 src/univariates.jl                       |  2 +-
 test/runtests.jl                         |  1 +
 test/weibull.jl                          | 27 ++++++++++++++++++++++++
 7 files changed, 75 insertions(+), 3 deletions(-)
 create mode 100644 docs/src/weibull.md
 create mode 100644 src/univariate/continuous/weibull.jl
 create mode 100644 test/weibull.jl

diff --git a/docs/make.jl b/docs/make.jl
index 658d477..22f0435 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -23,6 +23,7 @@ makedocs(;
         "Dependencies" => "set_optimize.md",
         "LogNormal" => "lognormal.md",
         "LogitNormal" => "logitnormal.md",
+        "Weibull" => "weibull.md",
         #"Details" => "z_autodocs.md",
     ],
 )
diff --git a/docs/src/weibull.md b/docs/src/weibull.md
new file mode 100644
index 0000000..bf73997
--- /dev/null
+++ b/docs/src/weibull.md
@@ -0,0 +1,18 @@
+```@meta
+CurrentModule = DistributionFits
+```
+
+# Weibull distribution
+
+The Weibull distribution has a scale and a shape parameter.
+It can be fitted using two quantiles.
+
+```jldoctest; output = false, setup = :(using DistributionFits,Optim)
+d = fit(Weibull, @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.
+
diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index 6053ee7..36ccd90 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -1,3 +1,5 @@
+
+
 function fit(::Type{Exponential}, m::AbstractMoments)
     # https://en.wikipedia.org/wiki/Exponential_distribution
     n_moments(m) >= 1 || error("Need mean to estimate exponential")
@@ -22,5 +24,3 @@ function fit_mode_quantile(::Type{Exponential}, mode::Real, qp::QuantilePoint)
     θ = -qp.q/log(1-qp.p)
     Exponential(θ)
 end
-
-
diff --git a/src/univariate/continuous/weibull.jl b/src/univariate/continuous/weibull.jl
new file mode 100644
index 0000000..14698b3
--- /dev/null
+++ b/src/univariate/continuous/weibull.jl
@@ -0,0 +1,25 @@
+# function fit(::Type{Exponential}, m::AbstractMoments)
+#     # https://en.wikipedia.org/wiki/Exponential_distribution
+#     n_moments(m) >= 1 || error("Need mean to estimate exponential")
+#     return Exponential(mean(m))
+# end
+
+function fit(::Type{Weibull}, lower::QuantilePoint, upper::QuantilePoint)
+    # https://www.johndcook.com/quantiles_parameters.pdf
+    gamma = (log(-log(1-upper.p)) - log(-log(1-lower.p)))/(log(upper.q) -log(lower.q))
+    beta = lower.q / (-log(1-lower.p))^(1/gamma)
+    Weibull(gamma, beta)
+end
+
+# function fit_mean_quantile(::Type{Exponential}, mean::Real, qp::QuantilePoint)
+#     # only fit to mean
+#     fit(Type{Exponential}, AbstractMoments(mean))
+# end
+
+# function fit_mode_quantile(::Type{Exponential}, mode::Real, qp::QuantilePoint)
+#     # ignore mode (its always at 0)
+#     θ = -qp.q/log(1-qp.p)
+#     Exponential(θ)
+# end
+
+
diff --git a/src/univariates.jl b/src/univariates.jl
index ec2689b..13e50d7 100644
--- a/src/univariates.jl
+++ b/src/univariates.jl
@@ -66,7 +66,7 @@ const continuous_distributions = [
     # "triweight",
     # "uniform",
     # "vonmises",
-    # "weibull"
+    "weibull"
 ]
 
 for dname in discrete_distributions
diff --git a/test/runtests.jl b/test/runtests.jl
index 3fbfb51..af381f5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,6 +17,7 @@ const tests = [
     "lognormal",
     "logitnormal",
     "exponential",
+    "weibull",
 ]
 #tests = ["logitnormal"]
 
diff --git a/test/weibull.jl b/test/weibull.jl
new file mode 100644
index 0000000..01a83fd
--- /dev/null
+++ b/test/weibull.jl
@@ -0,0 +1,27 @@
+
+# @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
+    D = Weibull(1,1)
+    #Plots.plot(D)
+    median(D)
+    qpl = @qp_m(median(D))
+    qpu = @qp_u(quantile(D, 0.95))
+    d = fit(Weibull, qpl, qpu);
+    @test d == D
+    #
+    qpl = @qp_m(0.5)
+    qpu = @qp_uu(5)
+    d = fit(Weibull, qpl, qpu);
+    @test median(d) == 0.5
+    @test quantile(d, 0.975) ≈ 5
+    #Plots.plot(d); 
+end;
+