From 928c74754b2e0f184e41f437905e887887956867 Mon Sep 17 00:00:00 2001
From: amyoshino <adriano.my@gmail.com>
Date: Wed, 28 Jun 2023 17:21:51 -0400
Subject: [PATCH] adding icdf functions for moyal, gumbel, triangular and
 weibull

---
 pymc/distributions/continuous.py       | 42 ++++++++++++++++++++++++++
 tests/distributions/test_continuous.py | 25 +++++++++++++++
 2 files changed, 67 insertions(+)

diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index fa4c056f93..3484a3b79a 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -2546,6 +2546,16 @@ def logp(value, alpha, beta):
             msg="alpha > 0, beta > 0",
         )
 
+    def icdf(value, alpha, beta):
+        res = beta * (-pt.log(1 - value)) ** (1 / alpha)
+        res = check_icdf_value(res, value)
+        return check_parameters(
+            res,
+            alpha > 0,
+            beta > 0,
+            msg="alpha > 0, beta > 0",
+        )
+
 
 class HalfStudentTRV(RandomVariable):
     name = "halfstudentt"
@@ -3089,6 +3099,20 @@ def logcdf(value, lower, c, upper):
             msg="lower <= c <= upper",
         )
 
+    def icdf(value, lower, c, upper):
+        res = pt.switch(
+            pt.lt(value, ((c - lower) / (upper - lower))),
+            lower + np.sqrt((upper - lower) * (c - lower) * value),
+            upper - np.sqrt((upper - lower) * (upper - c) * (1 - value)),
+        )
+        res = check_icdf_value(res, value)
+        return check_parameters(
+            res,
+            lower <= c,
+            c <= upper,
+            msg="lower <= c <= upper",
+        )
+
 
 @_default_transform.register(Triangular)
 def triangular_default_transform(op, rv):
@@ -3177,6 +3201,15 @@ def logcdf(value, mu, beta):
             msg="beta > 0",
         )
 
+    def icdf(value, mu, beta):
+        res = mu - beta * pt.log(-pt.log(value))
+        res = check_icdf_value(res, value)
+        return check_parameters(
+            res,
+            beta > 0,
+            msg="beta > 0",
+        )
+
 
 class RiceRV(RandomVariable):
     name = "rice"
@@ -3733,6 +3766,15 @@ def logcdf(value, mu, sigma):
             msg="sigma > 0",
         )
 
+    def icdf(value, mu, sigma):
+        res = sigma * -pt.log(2.0 * pt.erfcinv(value) ** 2) + mu
+        res = check_icdf_value(res, value)
+        return check_parameters(
+            res,
+            sigma > 0,
+            msg="sigma > 0",
+        )
+
 
 class PolyaGammaRV(RandomVariable):
     """Polya-Gamma random variable."""
diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py
index dccf6afeff..0e21e81d2f 100644
--- a/tests/distributions/test_continuous.py
+++ b/tests/distributions/test_continuous.py
@@ -207,6 +207,12 @@ def test_triangular(self):
             lambda value, c, lower, upper: st.triang.logcdf(value, c - lower, lower, upper - lower),
             skip_paramdomain_outside_edge_test=True,
         )
+        check_icdf(
+            pm.Triangular,
+            {"lower": -Rplusunif, "c": Runif, "upper": Rplusunif},
+            lambda q, c, lower, upper: st.triang.ppf(q, c - lower, lower, upper - lower),
+            skip_paramdomain_outside_edge_test=True,
+        )
 
         # Custom logp/logcdf check for values outside of domain
         valid_dist = pm.Triangular.dist(lower=0, upper=1, c=0.9, size=2)
@@ -704,6 +710,13 @@ def test_weibull_logcdf(self):
             lambda value, alpha, beta: st.exponweib.logcdf(value, 1, alpha, scale=beta),
         )
 
+    def test_weibull_icdf(self):
+        check_icdf(
+            pm.Weibull,
+            {"alpha": Rplusbig, "beta": Rplusbig},
+            lambda q, alpha, beta: st.exponweib.ppf(q, 1, alpha, scale=beta),
+        )
+
     def test_half_studentt(self):
         # this is only testing for nu=1 (halfcauchy)
         check_logp(
@@ -780,6 +793,11 @@ def test_gumbel(self):
             {"mu": R, "beta": Rplusbig},
             lambda value, mu, beta: st.gumbel_r.logcdf(value, loc=mu, scale=beta),
         )
+        check_icdf(
+            pm.Gumbel,
+            {"mu": R, "beta": Rplusbig},
+            lambda q, mu, beta: st.gumbel_r.ppf(q, loc=mu, scale=beta),
+        )
 
     def test_logistic(self):
         check_logp(
@@ -863,6 +881,13 @@ def test_moyal_logcdf(self):
         if pytensor.config.floatX == "float32":
             raise Exception("Flaky test: It passed this time, but XPASS is not allowed.")
 
+    def test_moyal_icdf(self):
+        check_icdf(
+            pm.Moyal,
+            {"mu": R, "sigma": Rplusbig},
+            lambda q, mu, sigma: floatX(st.moyal.ppf(q, mu, sigma)),
+        )
+
     def test_interpolated(self):
         for mu in R.vals:
             for sigma in Rplus.vals: